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I. INTRODUCTION 



Direction Finding/Polarization 
Estimation — Dipole and/or 
Loop Triad(s) 



KA1NAM THOMAS WONG, Senior Member, IEEE 
Chinese University of Hong Kong 



This paper shows 1) how measurement of the three Cartesian 
components of the electrical-field or magnetic-field suffices 
for multisource azimuth/elevation direction finding and 
polarization estimation, and 2) how the vector cross-product 
direction-of -arrival estimator is fully applicable even when the 
dipole triad is arbitrarily displaced from the loop triad. 
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A series of recent papers, [3-12, 14-18, 20, 21] 
among others, investigates the use of collocated 
six-component electromagnetic vector sensors for 
diversely polarized dkection-bf-arrival (DOA) 
estimation. A collocated six-component vector 
sensor consists of three identical and collocated 
but orthogonally oriented electrically short dipoles 
(called a dipole triad) plus three identical collocated 
but orthogonally oriented magnetically small loops 
(called a loop triad). All six component-antennas are 
spatially collocated in a point-like geometry. The 
collocated six-component vector sensor offers the 
following advantages: 1) the polarization diversity 
among the vector sensor* s component antennas allows 
that incident sources to be separated on account of 
their polarization differences in addition to their 
azimuth/elevation angular differences; 2) the spatial 
collocation of all component antennas in the vector 
sensor means no spatial phase delay in the vector 
sensor steering vector; hence, near-field sources 
may be located by an individual vector sensor as 
well as far-field sources; and 3) in a multisource 
scenario, each source's three Cartesian direction 
cosine estimates (and thus each source's azimuth 
angle estimate and the elevation angle estimate) are 
automatically paired without further post-processing. 
Theoretical performance bounds for direction finding 
using the collocated six-component vector sensors 
have been defined and derived in [3, 5], 

A variety of eigenstructure-based direction finding, 
polarization estimation and tracking schemes [6-12, 
14-18, 20, 21] have deployed these collocated 
six-component vector sensors in diverse array 
configurations for various signal scenarios using the 
vector cross-product DOA estimator. This vector 
cross-product DOA estimator exploits all six Cartesian 
components of the incident electromagnetic field 
to estimate the kth source's amplitude-normalized 
Poynting vector p A , which, in turn, gives estimates 
of the source's elevation angle B k (measured from the 
positive z-axis) and the azimuth angle </> k (measured 
from the positive x-axis): 



def 



LP*. 



11*11 11*111" 







"sinflfc 


cos^" 






sind k 


sin<^ 






.cos 6^ 





(1) 
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where * denotes complex conjugation, e k and h^, 
respectively, denote the fcth source's electric-field 
vector and magnetic-field vector, u k , v k , and 
w k> respectively, represent the fcth source's 
direction-cosines along the *-axis, the y-axis, and 
the z-axis. This vector cross-product DOA estimation 
approach complements the customary interferometry 
direction finding approach, which estimates the 
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spatial phase delay among the data sets collected at 
physically displaced antennas. 

The six-component electromagnetic vector-sensor, 
however, suffers from mutual coupling between its 
dipole triad and its loop triad. If only one triad (but 
not both) is deployed or if some significant distance 
separates the two triads, then the above mentioned 
coupling problem may be mitigated. Moreover 
deploying only one triad reduces antenna and receiver 
electronic hardware costs, simplifies algorithmic 
complexity, completely avoids the inter-dipole/loop 
triad mutual coupling problem, and eliminates the 
need to synchronize the phase between the dipole triad 
and the loop triad. 

This work shows 1) how a dipole triad by 
itself suffices for multisource azimuth/elevation 
direction finding and polarization estimation, 
2) how a loop triad by itself suffices for multisource 
azimuth/elevation direction finding and polarization 
estimation, and 3) how the vector cross-product DOA 
estimator remains fully applicable for a pair of dipole 
triad and loop triad spatially displaced by an arbitrary 
and (possibly) unknown distance (rather than being 
collocated). The above contrasts with the case of a 
collocated pair of perpendicularly oriented dipoles 
(commonly used in mobile communications), which 
is insufficient for exact estimation of the sources' 
arrival angles. In [19], however, it is shown that two 
horizontally oriented and magnetically small loops 
(plus an optional vertically oriented short dipole) 
can approximately estimate the sources' azimuth 
angles. Such an antenna array set-up matches the 
polarization of the wireless handset transmitter's 
strong vertical electric-field and can be used to 
retrofit dumb antenna-array receivers at the cellular 
base-station for downlink transmission beamforming. 

II. MATHEMATICAL MODELS OF STEERING 
VECTORS 

The fcth unit-power completely polarized 1 
transverse electromagnetic wave, having traveled 
through a homogeneous isotropic medium, is 
characterized by the electric-field vector e k and the 
magnetic-field vector h*, expressible in Cartesian 
coordinates as [3, 4]: 

^z( 9 k^Vk) 
hxWk^kylkyVk) 

h y( e k><t>k>yk>v k ) 
K^ik) 



CI* 



dcf 



COS<£ 4 COS0 A 

sin <p k cos $ k 
-sin*?* 
-sin<j> k 

COS(f> k 

0 



— sin0£ 
cos<f> k 
0 

-cos^cos0 A 



SU10* 



sin %e jrik "I 
cos 7* J 



dcf 
= gfc 



(2) 



where 0 < 9 k < n/2 denotes the signal's elevation 
angle measured from the vertical z-axis, 0 < <f> k < 2n 
symbolizes the azimuth angle measured from the 
positive x-axis, 0 < 7* < 7r/2 refers to the auxiliary 
polarization angle, and -7r < r) k < it represents the 
polarization phase difference. Note that 0(0*, <f> k ) 
depends only on the angular parameters, whereas g* 
depends only on the polarizational parameters. Also, 
Kll 55 \W\\ = 1 for all values of (O ki <t> k yif k ,Vk)- 

The dipole triad has the steering vector a fc = e*. 
The loop triad has the steering vector a* = h*. A 
dipole triad plus a displaced but identically oriented 
loop triad (with the loop triad located at (d x7 d y1 d z ) 
relative to the dipole triad) has the steering vector 



a = f ** I 



where q H (u k ,v k ) is defined as ^-W^-Mr****)/* 
and represents the spatial phase-factor relating the 
measurement at the loop triad to that at the dipole 
triad. A denotes the signals' wavelength. For reference, 



a, = 



LhJ 



for the collocated six-component vector sensor. 

The dipole triad, the loop triad, or the 
dipole-triad-plus-loop-triad displaced pair may be 
used as a multicomponent element in a multielement 
array for direction finding and polarization estimation, 
just as the collocated six-component vector sensor 
has been used in [3-12, 14-18, 20, 21]. Specifically, 
in [8] a single collocated six-component vector 
sensor (along with a temporai-invariance version 
of ESPRIT 2 ) to estimate the arrival angles and 
polarization states of multiple pure tones, while in 
[10, 21] the same is done for multiple frequency-hop 
(FH) signals. In [6, 7, 17, 18] a number of collocated 
six-component vector sensors are employed in a 
sparse (thinned) LxM rectangular array without 
incurring any cyclic ambiguity in the final estimates 
of the sources' Cartesian direction cosines by the 
use of a spatial-invariance version of ESPRIT. In 



'Partially polarized or unpolarized sources may be handled using 
the techniques presented in [12]. 



2 ESPRIT [2] is a closed-form eigenstructure-based parameter 
estimation technique that requires the data to possess certain 
"invariance" structures. 
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[9, 20], arbitrarily and irregularly spaced three- 
dimensional arrays of collocated six-component 
vector sensors are proposed in conjunction with a 
MUSIC 3 -based algorithm that needs no initial coarse 
estimates to start off MUSIC'S iterative search. In 
[11, 16], it is shown how the above arbitrarily spaced 
collocated six-component vector sensors are handled, 
using a polarization-invariance version of ESPRIT, 
when their locations are unknown. In all above 
mentioned schemes, the collocated six-component 
vector sensor may be substituted by the dipole triad, 
or the loop triad, or the dipole-triad-plus-loop-triad 
displaced pair; and this paper shows how. Without 
reciting the algorithmic details, all aforementioned 
schemes involve eigenstructure-based parameter 
estimation algorithmic steps that lead to an estimate 
of each incident source's steering vector, correct to 
within an unknown complex-value scalar c k . That 
is, available somewhere in each algorithm is the 
estimate i k « c k a k . (The approximation becomes 
a straight equality in noiseless or asymptotic 
cases.) The question is whether a* suffices to 
unambiguously estimate the &th source's arrival angles 
and polarization states when a* corresponds to the 
steering vector of a dipole triad, a loop triad, or a 
dipole-triad-plus-loop-triad displaced pair. The answer 
is "y es " and the following section shows how. 

III. ESTIMATION FORMULAS FOR THE ARRIVAL 
ANGLES AND POLARIZATION PARAMETERS 

The key observation is that ||e fc || = ||h k || = 
!IP*II = 1 for all (0*,0*'7*» T fc)' This means that the 
uncertainty in due to \c k \ may be removed by 
amplitude-normalizing a* m me cases of the dipole 
triad and the loop triad to produce a unity Frobenius 
norm. Algebraic and trigonometric manipulations on 
this normalized a*/||aj then produce five real- valued 
nonlinear equations, leading to unambiguous 
estimation of {0 k , <f> k , 7 fc , r} k }. 

A. For Dipole Triad 

If dipole triads are deployed, & k = c k e k under 
noiseless or asymptotic conditions. Hence, 

" - cos 0 k cos <f> k sin7 t + cos 7* sin <f> k cos tj 4 

- cos 6 k sin 4> k sm^ k — cos 7 t cos <f> k cos t) k 

sin0 t sin7 t 

" - cos 7 t sin <f> k sin r) k " 

+ j cos 7 t cos <f> k sin r) k (3) 

0 

where [•],■ refers to the ith element in the bracketed 
vector and L denotes the angle of the ensuing entity. 



The above array manifold model has not accounted 
for the dipoles' mutual coupling effects; however, 
good isolation and balance among the dipoles can 
minimize intratriad mutual coupling and renders 
the above array manifold model a very realistic 
approximation The above expression shows that 
the real and imaginary parts of e~ J ^** ]i & k /\\i k \], as 
five real-valued separate entities, producing for each 
incident source five separate nonlinear real- valued 
equations relating the four unknown signal parameters 
[B kf <t> k ,% 9 ri k }. Manipulation of these equations gives 



if lm{[SL k e-j^h] x } <o 



, m _ t -Im{[y-^ H} 

-ImUy-^b],} Im{[V -^ lj]l} > 0 
Im{[a t e-^l«»b) 2 } Ui J|/ ~ 

(4) 



-./II & \ 

I II Re{[a,] 1 e-> z l i »lJ}cos5 k + Re{[a t ) I «--' z Wj}sui^ J 

(5) 



* \ sin^ J 



(6) 
(7) 



Thus, the unknown complex-valued scalar ambiguity 
c k in a t needs to cause no ambiguity in the estimation 
of the arrival angles and polarization parameters. 



B. For Loop Triad 

If loop triads are deployed, SL k = c k h k under 
noiseless or asymptotic conditions. Hence, 



' - cos 9 k cos <j> k cos 7 t + sin 7 t sin <j> k cos r] k ' 
. - cos 9 k sin <f> k cos7 k - sin *y k cos <j> k cos rj k 
sin^cos7 t 



+ j 



-sin7 t sin0 t sin77 t 
sin7 t cos</> t sinT/ t 
0 



(8) 



3 MUSIC [1] is an iterative eigenstructure-based parameter 
estimation technique. 



The above array manifold model has not accounted 
for the loops' mutual coupling effects; however, 
good isolation and balance among the loops can 
minimize intratriad mutual coupling and renders 
the above array manifold model a very realistic 
approximation. The z-axis component of h k is always 
real in value regardless of the values of (0,0,7,77). 
With the left-hand-side of (8) already available, (8) 
produces for each incident source five nonlinear 
real-valued equations relating the four unknown 
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Fig. I. CRBs for various multicomponent antennas versus SNR. There are two incident narrowband completely polarized uncorrelated 
sources with (p v v v ti v y x ) = (0.41,0.3 l,7r/2,jr/4) and (k 2 »v 2 ,7?2,7 2 ) = (0.59,0.49, -7r/2,ir/4). 200 snapshots used at each SNR value. 



signal parameters. Manipulation of these equations 
gives 



where x denotes the vector cross product. Hence, 
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where I 3 symbolizes a 3 x 3 identity matrix, and 0 3 
refers to a 3 x 3 zero matrix. Hence, 



C For Displaced Dipole-Triad-Plus-Loop-Triad Pair 
Recalling that 

a * = C *L(«*.v t )hJ 

under noiseless or asymptotic conditions and that 
ll«H(«t." t )|| = 1. 

def - 



0 k = ai*sin(0pj? + [ftjf) = arccos([pJ 3 ) (16) 
k = arctan^Jj/IpJi) (17) 

(18) 



7. = arctan^—- 



where 



(20) 
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(21) 



IV COMPARATIVE CRAMER-RAO BOUND 
PERFORMANCE 

Fig. 1 plots the Cramer-Rao bound (CRB) 
[3] versus the signal-to-noise ratio (SNR) for the 
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several multicomponent antenna types discussed 
in this work under a scenario involving two 
closely, spaced, completely polarized, narrowband, 
uncorrelated far-field sources under additive 
white Gaussian noise. The source parameters are 
with (wi.Vpnj/Y!) = (0.41, 0.31, 7r/2,7r/4) and 
(« 2 ,v 2 ,7> 2 ,7 2 ) = (0.59,0.49,-7r/2,7r/4). There are 200 
snapshots used at each SNR value. 

The dipole triad and the loop triad has an exactly 
identical CRB curve, which is 3 dB higher than that 
of the collocated six-component vector sensor. This is 
expected because the former has half the number of 
component-antennas as the latter. The dash-dot and 
dash curves all refer to the dipole-triad-plus-lo op-triad 
displaced pair case, revealing the effects of the 
intertriad displacement axis' length and orientation. 
The two dash-dot curves and the top dash curve 
all correspond to an half-wavelength intertriad 
axis, whereas the bottom dash curve refers to a 
five-wavelength intertriad axis. The top dash-dot 
curve has the intertriad displacement aligned along 
the jc-axis; the lower dash-dot has the intertriad 
displacement aligned along the y-axis; and the 
upper dash curve has the inter-triad displacement 
45° from the x-y axes. It may be observed that a 
longer intertriad axis gives better estimates because 
of the larger geometric aperture that results. The axis 
orientation that optimizes estimation performance is 
the orientation onto which the incident sources project 
the farthest angular separation. In the present signal 
scenario, the 45° orientation (clockwise from the 
positive *-axis) gives the two sources a wider angular 
separation than when the intertriad axis is aligned 
either along the *-axis or along the y-axis. 
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Abstract: 

A novel eigenstructure scheme of extended-aperture spatial diversity and 
polarizational diversity is proposed for 2D arrival angle estimation. This 
innovational method involves uniformly spaced triads (or pairs) of electric 
dipoles or magnetic loops, spaced much farther apart than a half-wavelength: 
An ESPRIT-based step produces arrival angle estimates that suffer cyclic 
ambiguity due to the extended spacing. Then, a closed-form integer search is 
performed over a MUSIC null spectrum for all finite number of possible 
direction cosines in the cyclic ambiguity using the aforementioned cyclically 
ambiguous estimates. This approach produces highly accurate, yet 
unambiguous, azimuth and elevation angle estimates. This novel method 
facilitates (1) multipath de-correlation (independent fading), (2) coherent 
summation of multipaths from the same source, and (3) down-link 
beamforming from the base-station to the mobile 
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ABSTRACT 

A novel eigenstructure scheme of extended-aperture spatial 
diversity and polarizational diversity is herein proposed for 
2D arrival angle estimation. This innovational method in- 
volves uniformly spaced triads (or pairs) of electric dipoles 
or magnetic loops, spaced much farther apart than a half- 
wavelength. An ESPRIT-based step produces arrival angle 
estimates that suffer cyclic ambiguity due to the extended 
spacing. Then, a closed-form integer sear ch is performed 
over a MUSIC null spectrum for all finite number of possi- 
ble direction cosines in the cyclic ambiguity, using the afore- 
mentioned cyclically ambiguous estimates. This approach 
produces highly accurate, yet unambiguous, azimuth and 
elevation angle estimates. This novel method facilitates (1) 
multipath de-correlation (independent fading), (2) coher- 
ent summation of multipaths from the same source, and (3) 
down-link beamforming from the base-station to the mobile. 



1. INTRODUCTION 

Spatial diversity and polarizational diversity improve wire- 
less link performance. Spatial diversity and polarizational 
diversity are premised on the expectation that if one channel 
is in a deep fade, the multipaths may be constructive inter- 
fering in another channel. Note that field measurements J l] 
performed in a typical urban cellular environment have in- 
dicated that the horizontal and vertical polarizations have 
nearly equal power and are effectively uncorrelated. 

In frequency- selective slow-fading environments (slow- 
fading implies constant fading characteristics over a short 
observation period), spatial diversity and polarizational di- 
versity also allow various time-delayed multipaths to be re- 
solved based on their distinct polarizational states as well as 
their spatial arrival angles. These resolved multipaths may 
then be constructively summed to maximize the signal-to- 
noise ratio. Moreover, estimates of these sources' up-link 
arrival angles facilitate down-link beamforming from the 
base-station to the mobile. Spatial diversity is realized by 
deploying multiple spatially displaced antennas, typically at 
the base station. Polarization diversity also allows signals 
impinging from the same direction to be resolved on the 
basis of their different polarizational states. As the most 
significant scattcrers are typically close to the mobile, the 
base station antenna array aperture needs to extend at least 
over ten wavelengths to attain independent fading. 

The most straightforward way to realize a larger array 
aperture is to use more antennas, but at increased hardware 
and software costs. If it is desired to estimate the impinging 
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signals' arrival angles (for down-link beamforming and/or 
multipath constructive summation in frequency- selective 
fading as explained above) using the highly popular close- 
form computationally efficient eigenstructure method of ES- 
PRIT (Estimation of Signal Parameters via Rotational In- 
variance Techniques) [4], then spacing adjacent array ele- 
ments non-uniforraly over half a wavelength would generally 
violate ESPRIT's requirement of two identical but trans- 
lated subarrays. Uniform inter-element spacing exceeding 
half a wavelength, on the other hand, would lead to a set of 
ambiguous direction-cosine estimates. That is, the one-to- 
one correspondence between ESPRIT's eigenvalues and the 
direction-cosines would no longer exist when the uniform 
inter-element spacing exceeds half a wavelength. With no 
a-priori information, this ambiguity cannot be resolved us- 
ing unpolarized scalar sensors. 

This paper presents a direction-finding and polarization 
estimation scheme deploying a sparse array of either iden- 
tical triads or identical pairs of electric dipoles of magnetic 
loops to extend uniform inter-element spacing beyond a 
half- wavelength while resolving the aforementioned cyclic 
ambiguity. 

2. BASIC PRINCIPLES UNDERLYING NEW 
ALGORITHM 

Three co-located but orthogonally oriented electric dipoles 
would measure all three electrical-field components of the 
incident electromagnetic field separately. Similarly, three 
co-located but orthogonally oriented magnetic loops would 
measure all three magnetic-field components separately. To 
illustrate the basic principles of the method, the algorithm 
is described employing anLxM rectangular array of iden- 
tical triads or pairs of electric dipoles or magnetic loops. 
The inter- triad or inter-pair spacing A is assumed to be 
much greater than a half- wavelength to effect a large ar- 
ray aperture and correspondingly achieve highly accurate 
source direction estimates. One matrix-pencil is formed by 
considering the (L- I) x M triads or pairs on the first to the 
(L — i)-th rows as one sub-array and the (Ft — 1) X M triads 
or pairs on the second to the L-I.li rows as the other sub- 
array. This matrix-pencil has a spatial invariance along the 
x-axis and can yield estimates of the direction- cosines = 
sin*?* cos^jt, k = 1,..., A'}, where 0 < 6k < tt/2 is the 
signal's elevation angle measured from the vertical Zraxis, 
0 < 4>k < 2tt is the azimuth angle. The second matrix- pencil 
isTormed by considering the (M -t)xL triads or pairs on 
the first to the (M — l)-th columns as one sub-array and the 
the (M — 1) x L triads or pairs on the second to the M-th 
rows as the other sub- array. This matrix-pencil has a spatial 
invariance along the y-axis and can yield estimates of the 
direction- cosines {vj — sin0> sin^j, j — I, .... A'}. How- 
ever, when the inter- triad or inter- pair spacing is greater 
than a half- wavelength there is an integer multiple of 2jt 
ambiguity in the phase of each eigenvalue generated in the 
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final stage of TLS-ESPR1T. This leads to ambiguous DOA 
estimates equi-spaced by A/ A in the interval — 1 < Uk < 1 
(or — 1 < Vk < 1 as the case may be.) If one could re- 
solve this ambiguity, the resulting DOA estimate would be 
highly accurate due to the large aperture. The key idea of 
this paper is to determine the minimum of the MUSIC null 
spectrum by stepping through all finite number of possi- 
ble direction-cosine values, using the cyclically ambiguous 
estimates. 

3. SPARSE ARRAY DATA MODEL 

The propagation model is transverse electromagnetic plane- 
waves traveling through a non- conductive homogeneous 
isotropic medium. Under these conditions, the wavefront 
incident upon the array has an electric- field vector, e, and 
corresponding magnetic-field vector, h, that can be ex- 
pressed in Cartesian coordinates in matrix form as: [2] 

sinik cosOk cosfa e JT? * — cosyk sin<f>k 
sin-fk cosOk siruj>k e J71k -f cos 7* cos<f>k 
def -5*n7fc sinOk e^ k 

—cosjk cosOk cos<f>k — sinyk sinfa e JT? * 
-cosfk CosOk sinfa + sinyk cos<t>k e )Vk 
_ cosfksinOk 



noise vector at the (J, m)-th triad or pair, © is the Kronecker 
product, A is the wavelength associated with the center fre- 
quency, uf Ci of the passband of the front end bandpass filter. 

Furthermore, a* d = a(0fc, <f>u, 7*, f?jk ) is the 7x1 steering vec- 
tor for the fc-th source sensed at each triad or pair. That 
is, for triads of dipoles: 



a fc = 



(1) 



COS ^Jfc cos Ok 

sin 4>k cos 0k 

— sin $k 

— sin <f>k 
cos 4> k 

0 



— sin <f>k 
cos <f>k 

0 

• COS <f>k cos $k 

- sin <f>k cos 0* 
sin Bk 



sin7fe e™* 
cos 7* 



(2) 



where 0 < 7 < tt/2 is the auxiliary polarization angle, — n < 
» < 7T is the polarization phase difference, 0 < 6 < ir is the 
signal's elevation angle measured from the vertical a-axis, 
0 < <t> < 2tt is the azimuth angle, Z 0 is the transmission 
medium's intrinsic impedance. 

A* co-channel signal sources are assumed to be in the far 
field and narrowband in the sense that the bandwidth is 
very small relative to the carrier frequency. Given a uni- 
form rectangular array of the aforementioned identical tri- 
ads or identical pairs of electric dipoles or magnetic loops, 
the JLM x 1 snapshot vector (where J = 3 for an array of 
triads and J — 2 for an array of pairs) at time t may be 
expressed as 



(3) 



Z M = **(0<b(«0 ® <ly( v fc) ® a (#fci ^fci 7k» Ik) + n{t) 

= As(i) + n(t) 
where A is the JLM X K matrix: 

A d = [q*(ui ) ® <\ y (v\ ) <g> a(0i , ^i, 71 , 171 ), • • ■ , 

<\x(u K ) <8> <\y{vi<) ® a(0/c, <t>K, 7/c, r) K )} (4) 

where 









' ni,,(t) " 


/.v def 




; n(t) = 
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q*(«*) = 



j2wlt\*u k 



e 



e "X — 



e * 



s*(t) is the complex baseband signal comprising the k-th 
signal arrival, rti, m (t) is the J x 1 additive zero- mean white 



e Vk 



cos <t>k cos 0* - sin <j>k 
sin <^fc cos Ok cos 
- sin Bk 0 



T sin 7* e JVk 1 
L cos 7* J 



(5 



For triads of magnetic loops: 



afc = 



— sin 4>k — cos <j>k cos Bk 
cos 0* — sin <f>k cos 0* 
0 sinflfc 



[ sin 7* e" k ] 
L cos 7* J 



Furthermore, if dipole pairs or loop pairs are used instead 
of dipole triads or loop triads, then a* becomes 2 x 1, with 
the appropriate Cartesian component deleted from (5) or 
(6). It would also be possible to have co-located mix of 
dipoles and loops. For example, a dipole oriented along the 
ar-axis, another dipole along the y-axis plus a loop along 
the y-axis. When all six electromagnetic- field components 
are separately measured, then a simpler procedure alter- 
native to the present algorithm may be used [5]. Note 
that in any case, 0* depends only on the sources 1 spa- 
tial angular location and g* depends only on the incident 
signals' polarization states. The maximum number of un- 
correlated sources that the present algorithm can handle 
with such a L x M rectangular array of vector-sensors is 
min{J(L - l)M - 1, JL(M - 1) - 1}. 

With a total of iV > K snapshots taken at the dis- 
tinct times {i n , n = 1, . . . , AT}, the present direction- finding 
problem is to determine {0k^4 > k, k = from the 

JLM x N data set: 



r Zi.i 

Zl.Af 



Zi,i(£;\r) 



zi,m(*i) * ■ zi,m(£.w) 

zt.i(ii) ... zl,i(£n) 
zl,m(*i) ... zl,m(^) 



(7) 



where each of the LM sub- matrices zi, m of size J X N cor- 
responds to data measured at the m)-th vector-sensor. 

4. 2D ANGLE ESTIMATION WITH A SPARSE 
ARRAY OF ELECTRIC DIPOLES OR A 
SPARSE ARRAY OF MAGNETIC LOOPS 

The first step in ESPRIT is to compute the K (JLM x 
l) signal-subspacc* eigenvectors by eigen-decomposing the 
JLM x JLM data covariancc matrix ZZ T . Although two 
distinct matrix- pencil pairs are to be formed, only one 
eigen-decomposition of this data eovariance matrix is neces- 
sary because both matrix-pencils are formed from the same 
JLM x .V data set. Let E* be the JLM x K matrix com- 
posed of the K eigenvectors corresponding to the A"* largest 
eigenvalues of the JLM x JLM sample covariancc matrix 
Kzz = ZZ W . Due to the invariance associated with the two 
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subarrays along the x- axis, it follows from the snapshot 
model in (3) that Es asymptotically approaches Es =s AT, 
where T is an unknown K x K non-singular matrix to be 
determined. T is necessarily non-singular because both A 
and Es are full- rank matrices. 

4.1. Deriving the low- variance but ambiguous es- 
timates of Ukl 
For the matrix pencil with spatial in variance along the x- 
axis, define E? as the first J(L — l)M rows of Es and E£ 
as the last J(L — \)M rows of Es, that is: 



3?= [o 



J(L-l)M 



J Es = a; j ^ 

-l)A#j Es = • 



(8) 
(9) 



where 



EJ= [Qjxf : Ij(L-i)A/j Es = AJT 
A 2 U = A, = [q<?>(ui ) © q*(f i ) ® , , 7 i , *h ), ■ - , 



(I)/ \ del 



e 



qt 2) («u) = 



and I, is a » x t identity matrix and O, is an a x i matrix 
whose elements are all zeroes. In noiseless or asymptotic 
cases, there exists a /v x K non-singular matrix *fc tt relating 
the two J(L - \)M x A* full-rank matrices E? and EJ [4]: 

E?* 11 = ES AiT u V = A?T U (10) 

=> * u = {(En H EJ t }~ , {(Ei 4 ) H E2} (11) 
=> tf" = (T u )- l * u T u (12) 

where <fr u is a AT x A' diagonal matrix whose diagonal ele- 
ments {[$ u ]fct = e j2 T A * Uk , fc =s 1, . . . , A'} constitute the 
eigenvalues of 4? tt . Because A* > ™ and -1 < u* < 1, 
there exists a set of cyclically related candidates for the 
estimation of u&: 

[^(-1-H fc )| <n m < [^(l-itt)] (13) 

where \x] is the smallest integer not less than x, \x\ is the 
largest integer not greater than and L refers to the princi- 
pal argument. Furthermore, the eigenvector corresponding 
to the eigenvalue [& v ]kk is the fc-th column of (T w )~ l , where 
T u — P 11 T and P u is an unknown permutation matrix. In 
more realistic cases where noise is present and when only 
a finite number of snapshots are available, all of the above 
estimates become only approximate. 

4.2. Deriving the low-variance but ambiguous es- 
timates of viti 
Similarly, for the matrix pencil with spatial in variance along 
the y-axis, define EJ and EJ as: 



Er = 



(14) 



Ej = 



h J 



e s = a;t 



(15) 



L-block-diagonal 



where 



A? = [q*(iii)®q£ 1) (vi)$a(0i,4ii7iitt).'-'. 
q*(u/c) <S> q v l) (v K ) ® a(0/c, <t>K> 7k, m)) 
Al = A\ & = [q a (ui ) ® q<, 2) (t/, ) 0 a(0, , , 71 , m ), ■ • ■ , 
q*(u*r) 0 q£ 2) (u jo) ® a(^A- t 4>k< 7a* , »?*•)] 



qy («*) = 



e * 



L e * 





r J 3 » A » a "* n 




e * 


(2) / \ def 






L e * J 



/,-block-diagonal 



fi - [l3 ( w-o : 0 3 h = [o 3 ! I 3(M _„] (16) 

In noiseless or asymptotic cases, the solution to the matrix 
equation E\* v = E? is * v = T- l * v T, where T v = P V T 
and P w is another unknown permutation matrix. Let the 
eigenvalues of ^ v be denoted i/>, j = 1,...,A\ It follows 
that 

^-[•Tfl j = l,... t A' (17) 

^ w = diag |e , '^ Ay1 ' 1 , e ,: ^ A * Va , • • ■ , | (18) 

and the corresponding eigenvectors are the columns of T" 1 . 
With A y greater than a half wavelength, the direction- 
cosine candidates of the j-th source relative to the t/-axis 

\^(-l-^)]<n v <[^(l-^)\ (19) 

(19) represents a set of of low- variance but ambiguous esti- 
mates of the direction-cosine relative to the y-axis. 

4.2.1. Pairing the Direction- Cosine Estimates: 
Note that different indices are used to enumerate the 
eigenvalues of and (i.e., the low-variance but cycli- 
cally ambiguous estimates of the direction-cosines k = 
1,...,A'} and j = l t ...,A'}) This is because even 
though 4i u and 3f have common eigenvector (in the noise- 
less or asymptotic cases), the ordering of the eigenvectors of 
ty v will generally be permuted relative to the ordering of the 
eigenvectors of ^ u due to the presence of the unknown per- 
mutation matrices P u and P . In other words, the same 
set of eigenvectors is ordered differently in (T 11 )" 1 as in 
(T*)" 1 . No mismatch, however, exists between /jjt and its 
corresponding eigenvector, which is simply the fc-th column 
of (T u ) -1 . Likewise, no mismatch exists between vj and its 
corresponding eigenvector, which is simply the j-th column 
of (T w ) _1 . Thus, fib may be paired with i/j from the same 
source by matching the orthogonal rows of T" with those 
of T w as follows. Let (jjt, k) denote the row- index of the 
matrix element with the largest absolute value in the fc-th 
column of the K x K matrix T'YT*)- 1 . Then n* and u 2k 
belong to the same source. Note tnat this pairing procedure 
involves minimum computation and requires no exhaustive 
searches. 
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4.3. MUSIC-Based Disambiguation 
After the overall element- space is decomposed into a K- 
dimensional signal subspace and a (LM — /tridimensional 
noise subspace, MUSIC (MUltiple Signal Classification) [2] 
constructs a null spectrum parameterized by the signal pa- 
rameters. Iterative search over this null spectrum for its 
K deepest nulls would produce estimates of the K sources* 
parameters. 

MUSIC's null spectium is constructed based on the or- 
thogonality between the signal subspace and the noise sub- 
space of the sample correlation matrix R**. That is, 
JL E n - Because E* and A span the same column space, 
A 1 E n - Recalling that u and v are both functions of 



arg min 



I En (<M«) ® q»(f) ® <t>< 7,.»7) 



where the cost function above may be rewritten as: 

= M«) ® q B («) ® (*(* , *)g(7. l))) H EnE^ 
[<b(«) 9 qy(w) ® (*(*. #)g(7, »?))] 
= [(q»(«) ® «b(»))* ® (S(7, r,) H *(0, </.)")] E„E* 
[q»(«) ® q B (f ) ® (*(*. 0)g(7- »/))] 
= g(7,V ? ) w M(tf,^)g(7.D) 

where 

M(tf, =' [(q,(«) ® q s (tO)" ® *(«, 

E„E? (q,(u) ® q,(u) ® $(0, <t>)] 



(20) 



The identity (M 1 M 2 )©(M3M < ) = (Mi ® M 3 )(M2® M 4 ) 
is used in the preceding step. Noting that M(0, <£) is a 2 x 2 
non-negative definite matrix, 



{ft 



7 . arg mm 



{smallest eigenvalue of M(0,<£)} (21) 



Thus, the direction- finding problem is again decoupled from 
the polarization estimation problem. 

Note that this MUSIC-based disambiguation step renders 
it unnecessary to use the vector cross-product estimator to 
derive any coarse reference estimates. Thus, it becomes 
no longer necessary to have all available all six electromag- 
netic components of the impinging wavefront. It would be 
sufficient to have available as few as only any two of the 
six electromagnetic components, because the resulting ar- 
ray manifold would still retain its one-to-one relation with 
the impinging source's direction cosine. The cyclically am- 
biguous set of direction cosines along the x-axis may be 
aligned with the set along the y-axis by methods such as 
that presented in [8]. 

While the present MUSIC step may be applied to arrays 
of any geometry, the uniformly spaced configuration used 
in the present algorithm minimizes the support region of 
MUSIC's iterative search from a continuous region {0 < 
$ < 7T, 0 < <p < 2ir} to a small finite set of discrete points: 

**("*). for f^(-l - Wk)] < «« < [^(1 - Wk)] . 

[^ M -, fc )l<n.<[^(l-.,j} 

5. SIMULATIONS 

Simulation results in Figures 1 to 4 verify the efficacy of the 
proposed extended aperture concept using sparsely spaced 
dipoles and/or loops. The wireless mobile fading channel 



is modeled as slow-fading and frequency-selective. For il- 
lustration, present simulations use a two-multipath model. 
The source is a BPSK sequence modulated by the raised- 
cosine waveform with f3 = 0.35. The two multipaths, with 
timing offset equals 0.3 of a symbol period, impinge upon a 
4x4 rectangular array of identically oriented dipole- triads 
equi-spaced at 8 wavelengths. Each dipole triad consists of 
three co-located but orthogonally oriented identical dipoles. 
The multipath parameters were: ui = 0.61, vi = 0.79, 71 = 
45°,tn = -90°, "Pi = 1, tii = 0.69,v 2 = 0.71,72 = 45°,^ = 
90°, V2 = 0.7079. 50 symbols are sampled at 10 samples per 
symbol in each experiment; and there are 100 experiments 
per data point. Spatial smoothing performed by forming 
four 3x3 dipole-triad subarrays. For comparison, the same 
scenario is used with a 7 x 7 half- wavelength spaced omni- 
directional scalar-sensor array, where spatial-smoothing [3] 
is performed using four 6x6 subarrays. A scalar sensor 
refers to the situation where each dipole-triad is replaced 
by a scalar- sensor which measures the Frobenius norm of 
the Poynting vector. This 7x7 half- wavelength spaced 
scalar-sensor array has an array manifold size comparable 
to the 4x4 ex tended- aperture dipole-triad array. 

Referring to figures 1 and 2, the present extended- 
aperture algorithm clearly offers significantly lower es- 
timation standard deviations and biases than the half- 
wavelength- spaced scalar-sensor array at all SNR. Because 
the two multipaths' direction- cosines are separated by 0.08, 
the two multipaths may be considered as resolved when 
their estimation standard deviations fall below 0.01. Using 
this criteria, the proposed scheme has a resolution thresh- 
old below -10 dB compared to about 0 dB for the half- 
wavelength spaced scalar-sensor array. In the SNR > 0 dB 
range where the standard deviation and bias of both arrays 
fall off rather linearly with increasing SNR, the proposed 
method offers almost one and a half orders of magnitude 
reduction in standard deviation and about an order of mag- 
nitude reduction in bias. 

In figures 3 and 4, the SNR becomes set at 0 dB while 
the muJtipaths' relative timing offset is varied. Again, the 
present extended-aperture algorithm clearly offers signifi- 
cantly lower estimation standard deviations and biases than 
the half- wavelength-spaced scalar- sensor array at all timing 
offsets. The scalar-sensor array, despite spatial smoothing, 
fails to resolve the two multi-paths for timing offset below 
0.3 of a symbol period. In contrast, the proposed scheme 
offers near identical error statistics to all timing offset. In 
the timing offset range greater than 0.5 of a symbol period, 
where the half-wavelength spaced scalar-sensor array's er- 
ror statistics level off, the proposed scheme outperform the 
half- wavelength spaced seal ar- sensor array, offering about 
an onler of magnitude reduction in standard deviation and 
bias. 

When all six electromagnetic field components are to be 
separately measured, (for example by a "vector sensor" 
comprising three orthogonally oriented dipoles plus three 
orthogonally oriented loops, all co- located in space,) then 
an alternative disambiguation approach is realizable oy ap- 
plying a vector cross-product between the electric-field com- 
ponents and the magnetic-field components. For details, 
please refer to [5,6]. An underwater acoustic analog of the 

t>resent algorithm is introduced by the present authors in 
7]. Another extended-aperture algorithm, applicable to 
any propagation medium, is presented also by the present 
authors in [8] using dual-size spatial invarianccs. 
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Figure 1: RMS standard deviation of {ur, 0i, «2, M, vs. 
SNR: 

10* 




o-HafS-Wavelength Spaced Scalar Seniors 
x-Spareefy Spaced DlpcJe Tralds 

ul-0.61, u2-0. 79, v 1-0.69, v2.a71.SpacinQ.16 haU-wavetengtns 



Figure 2: RMS Bias of {«i , Oi , u 2 , $2}, vs. SNR. All param- 
eters are identical with those in Figure 1. 
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Figure ,3: RMS standard deviation of {«i,t5i, u 2 , O2} vs. 
timing ofFset, SNR with respect to first multi-path is 0 dB, 
all other parameters same as in figure 1. 
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Figure 4: RMS Bias of ^u\,vi,u 2 , O2} vs. timing offset, all 
parameters are same as in figure 3. 
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Abstract: 

A novel ESPRIT-based 2D angle estimation scheme is proposed involving the 
use of a right-triangular array of three vector sensors spaced much farther 
apart than a half-wavelength. A vector sensor is composed of six co-located 
antennas distinctly measuring all six electromagnetic field components of a 
multi-component incident wavefield. Information on each source's respective 
electromagnetic field components is obtained by decoupling the signal 
eigenvectors via lower dimensional eigenvectors derived from TLS-ESPRIT. 
This facilitates estimation of each source's respective Poynting vector thereby 
enabling one to resolve the phase ambiguities in ESPRIT's eigenvalues when 
the intervector-sensor; spacing is greater than a half-wavelength. Simulations 
are presented showing the sample variance of the direction cosine estimates 
decreasing linearly on a log-log scale as the intervector-sensor spacing is 
increased from a half-wavelength to 30 wavelengths, with a factor of 80 
reduction in the latter case relative to the former case. The proposed scheme 
and attendant vector sensor array also outperform a uniformly-spaced array 
of scalar sensors with the same aperture and same number of component 
antennas whenever the intervector sensor spacing in the former case is 
greater than 3 half-wavelengths 



Index Terms: 

array signal processing direction-of-arrival estimation eigenvalues and eigenfunctions 
signal resolution ESPRIT based 2D angle estimation Poynting vector TLS-ESPRIT 
array resolution colocated antennas direction cosine estimates eigenvectors 
electromagnetic field components extended aperture vector sensor arrays half 
wavelength high accuracy 2D angle estimation intervector-sensor spacing 
multicomponent incident wavefield phase ambiguities right triangular array sample 
variance scalar sensors signal eigenvectors simulations uniformly spaced array 
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ABSTRACT 

A novel ESPRIT-based 2D angle estimation scheme is 
proposed involving the use of a right- triangular array of 
three vector sensors spaced much farther apart than a half- 
wavelength. A vector sensor is composed of six co-located 
antennas distinctly measuring all six electromagnetic Held 
components of a multi-component incident wavefield. In- 
formation on each source's respective electromagnetic field 
components is obtained bv decoupling the signal eigenvec- 
tors via lower dimensional eigenvectors derived from TLS- 
ESPRIT. This facilitates estimation of each source's respec- 
tive Poynting vector thereby enabling one to resolve the 
phase ambiguities in ESPRlTs eigenvalues when the inter- 
vector-sensor spacing is greater than a half- wavelength. 
Simulations are presented showing the sample variance of 
the direction cosine estimates decreasing linearly on a log- 
log scale as the inter- vector-sensor spacing is increased from 
a half-wavelength to 30 wavelengths, with a factor of 80 re- 
duction in the latter case relative to the former case. The 
proposed scheme and attendant vector sensor array also 
outperform a uniformly-spaced array of scalar sensors with 
the same aperture and same number of component anten- 
nas whenever the inter-vector sensor spacing in the former 
case is greater than 3 half-wavelengths. 



1. INTRODUCTION 

Array resolution and accuracy are limited by the effective 
size of the array aperture. In general, a larger array aper- 
ture yields correspondingly more accurate DOA estimates 
[1]. The array aperture may be enlarged by adding more 
antenna elements in the case of uniform half-wavelength 
spacing, by spacing the elements nonuniformly over a larger 
aperture, or by increasing the separation between elements 
in the case of uniform spacing. Adding more antennas 
has the obvious drawbacks of increasing hardware costs 
plus expanding the already considerable computational load 
required by eigenstructure methods. Nonuniform inter- 
element spacing would generally violate ESPRIT's require- 
ment of two identical but translated subarrays. The spatial 
version of the Nyquist Sampling Theorem also poses an up- 
per limit on the distance between elements in the case of 
uniform spacing without causing aliasing. Two identical 
sensors spaced over a half- wavelength apart will lead to a 
set of ambiguous direction-cosine estimates. With no a- 
priori information, this ambiguity cannot be resolved using 
unpolarized scalar sensors. 

A novel scheme is proposed employing a sparse array 

tThis research was supported by the Air Force Office of Sci- 
entific Research under grant no. F4 0620-95- 1-0367, the National 
Science Foundation under under grant no. MIPS- 9320890, and 
the Army Research Office's Focused Research Initiative under 
grant number DAAH04-95- 1-0246. 



of six-component electromagnetic vector-sensors [2] that is 
able to resolve the aforementioned ambiguity despite inter- 
vector sensor spacings much greater than a half- wavelength. 
Exceptional estimation accuracy is facilitated by a large ef- 
fective aperture and exploitation of differences in sources' 
respective polarizations. 

1.1. Basic Principles Underlying New Algorithm 

A vector-sensor consists of six spatially co-located antennas 
measuring all three electrical-field components and all three 
magnetic-field components of the incident wave-field [2]. 
The proposed array geometry is that of three vector-sensors 
located at (0,0), (A,0), and (0, A) in (x,y) Cartesian co- 
ordinates. The inter-vector sensor spacing A is assumed 
to be much greater than a half- wavelength to effect a large 
array aperture and correspondingly achieve highly accurate 
source direction estimates. A pair of identical vector-sensors 
effectively represent two identical six-element subarrays dis- 
placed in space thereby facilitating the use of ESPRIT [3]. 
TLS-ESPRIT is applied to the data from the vector sensors 
located at (0,0) and (A,0) to ultimately yield estimates of 
the direction cosine of each source relative to the x-axis, 

{u k — cos 4> k sin 9k , k = 1 /C}, where <j> k and 0 k are 

the azimuth and elevation arrival angles of the fc-th incident 
signal. Simultaneously, TLS-ESPRIT is applied to the data 
from the vector sensors located at (0, 0) and (0, A) to ulti- 
mately yield estimates of the direction cosine of each source 
relative to the j^axis, {v k = cos<f> k sin & k} k - 1, . . . , K}. 

However, when the inter-vector sensor spacing is greater 
than a half- wavelength there is an integer multiple of 2k 
ambiguity in the phase of each eigenvalue generated in the 
final stage of TLS-ESPRIT. This leads to ambiguous DOA 
estimates equi-spaced by A/ A in the interval — 1 < u k < 1 
(or -1 < v k < 1 as the case may be.) If one could resolve 
this ambiguity, the resulting DOA estimate would be highly 
accurate due to the large aperture. The key idea of this 
paper is to use the Poynting vector information for each 
source embedded in the output of each vector sensor to 
resolve this ambiguity. The lower dimensional eigenvectors 

§enerated in the final stage of TLS-ESPRIT are used to 
ecouple the element space signal eigenvectors to yield both 
the electric field vector and the magnetic field vector for 
each source. The cross product between these two vectors 
for a given source yields the Poynting vector for that source. 
Normalizing the Poynting vector for a given source to have 
unit length, its components are the direction cosines of that 
source relative to the z, y, and z axes. That is, 
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where u>* =s y/l — uj — v\, * denotes conjugation, and 7* 
and r?* are polarization parameters defined in Section 2. 
The idea is to use these as "coarse 11 estimates to resolve the 
integer multiple of XJA ambiguity in the DO A estimates 
yielded by TLS-ESPRIT's eigenvalues. 

The direction cosine estimates produced by the Poynt- 
ing vector estimation procedure are characterized as high 
variance but unambiguous in comparison to the direction 
cosine estimates extracted from the TLS-ESPRIT eigenval- 
ues which are characterized as low variance but ambiguous. 
Intuitively, the lower variance of the latter relative to the 
former is due to the fact that the direction cosine estimates 
extracted from the TLS-ESPRIT eigenvalues depend on the 
size of the array aperture as measured by A: the larger A, 
the smaller the variance of Z/ik/2ffA, where p k is a TLS- 
ESPRIT eigenvalue associated with the fc-th source. This 
is confirmed by simulations presented in Section 4. The 
relative high variance of the direction cosine estimates ex- 
tracted from the Poynting vector estimates is intuitively 
due to the fact that they are inherently extracted from in- 
formation provided by a single vector sensor which has no 
effective aperture. 

2. VECTOR SENSOR ARRAY DATA MODEL 

The propagation model is transverse electromagnetic 
planewaves traveling through a non-conductive homoge- 
neous isotropic medium. Under these conditions, the wave- 
front incident upon the array has an electric-field vector, 
e, and corresponding magnetic-field vector, h, that can be 
expressed in Cartesian coordinates as [2] 

e = e,v x + e v v y + e 2 v\ (2) 
= {siny cos9 cos<f> e JTf — cos 7 sin<i>)v x 
-f (stn7 cos9 sin<f> e 3 * 1 + cosy co$4>)$ v 
-siny sinB e iv ^ z 
h = M*+Mv+^^ 

= — (co*7 cosQ costf) + siny sin<j> e }v ) Z 0 ^ x 
—(cosy cosB sin<}> — sin 7 cos<j> e JT *) Z 0 v y 
+cosy sxn9 Z Q v, (5) 

where 0 < 7 < tt/2 is the auxiliary polarization angle, — tt < 
n < v is the polarization phase difference, 0 < 9 < ir is the 
signal's elevation angle measured from the vertical z-axis, 
0 < <f> < 2tt is the azimuth angle, Zo is the transmission 
medium's intrinsic impedance, and v is a unit- vector along 
the coordinate specified by its subscript. 

K co-channel signal sources are assumed to be in the far 
field and narrowband in the sense that the bandwidth is 
very small relative to the carrier frequency. Given a right- 
triangular array of three vector sensors with the geometry 
described previously, the 18 x 1 snapshot vector at time t 
may be expressed as 



z(t). 



:£>(t)e"" 

kml 



v(*,**)' 



a fc = 



(3) 
(4) 



'v(^,^)®a(^ t ^,7*.»7*)+n(t). (6) 



The first, middle, and last 6x1 sub-vectors of z(t) are 
the outputs of the six component-antennas comprising the 
vector sensor located at the Cartesian coordinates (0,0), 
(A,0), and (0, A), respectively, relative to the x — y plane. 
The various quantities are defined as follows. v($k<4>k) is 
the 3x1 manifold describing the relative phases between 
vector sensors for the k-th source due to the propagation 
delays 

1 



(7) 



where A is the wavelength associated with the center fre- 
quency, u> c , of the passband of the front end bandpass fil- 
ter. <g> is the Kronecker product, a* =* &(9k,4>k, 7*, is 
the 6x1 composite electromagnetic field vector for the &-th 
source sensed at each vector sensor, i.e., 

sin 7jt cosOk cos<f>k e 3r,k - cosyn sinfa 
sinykCos9ksin<f>ke' llk + cosykcos<t>h 
—sinykS\n9kC }rtk 

(-cosykC089hCOS<f>k - sinyksin<f>ke 3t)k )Z 0 
(-cosykCOs9 k sin<l>k + 8inykCos4> k e 3Vik )Z 0 



Sk(t) is the complex baseband signal comprising the fc-th 
signal arrival. Finally, n(t) is an 18 x 1 additive noise vector. 

3. 2D ANGLE ESTIMATION VIA AN 
EXTENDED APERTURE VECTOR 
SENSOR ARRAY 

3.1. Adapting ESPRIT to Vector Sensor Arrays 
ESPRIT [l] is applicable to the right triangular array of 
three vector sensors because any two vector sensors repre- 
sents a pair of identical six element subarrays. The first step 
in ESPRIT is to compute the signal eigenvectors. Let Es be 
an 18 X K matrix composed of the K "largest" eigenvectors 
of the 18 x 18 sample covariance matrix 



N 

•=1 



Due to the invariance associated with the two vector sensors 
along the y-axis as well as the invariance associated with 
the two vector sensors along the x-axis, it follows from the 
snapshot model in (6), that E5 asymptotically approaches 
the following form 



Es : 



Ei 




Ax 


E 2 


— AT = 


A x * u 


E 3 




Aiw" 



T, 



(9) 



where A 1 is the 6 x K matrix 

a, = [.„«„..., , £ ]. do) 

$ u and & v are the K x K diagonal matrices, 

^ = diag{e^ Aul ,e^ AU3 e*'* A **}, (11) 



= diag 



(12) 



respectively, and T is an unknown K x K matrix. 

It is easily shown that in the no noise or asymptotic cases 
the solution to the matrix equation fei'l FU = is ty u = 
T" I ^ U T. Let the eigenvalues of flf u be denoted /i*, k = 
1,...,K. It follows that /i fc = Vtk = e J "¥ A "\ Jb = 1, 
and the corresponding eigenvectors are the columns of T -1 . 
If the inter- vector sensor spacing A is greater than a half 
wavelength, the candidates for the direction cosine of the 
fc-th source relative to the x-axis are 

( i3 > 

where Z/U is the principal argument of fik in the range 
— n < l\ik < tt, fx] is the closest integer to x that is greater 
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than x, and |xj is the closest integer to x that is less than 
x. (13) represents a set of of low variance but ambiguous 
estimates of the direction cosine of the k-th source relative 
to the x-axis. 

Similarly, in the asymptotic or no noise cases the solution 

to the matrix equation Ei4f v = E 3 is * w = T~ l *"T. 
Let the eigenvalues of * p be denoted u it I as 1, .... K. It 

foUows that vt = = / = 1, and the 

corresponding eigenvectors are the columns of T _1 . With 
A greater than a half wavelength, the candidates for the 
direction cosine of the l-th source relative to the j/-axis are 

.(») A A f A Lv k \ ^ I A lu k \ 

(14) 

(14) represents a set of of low variance but ambiguous esti- 
mates of the direction cosine of the fc-th source relative to 
the {f-axis. 

Note that different indices are used to enumerate the re- 
spective eigenvalues of * tt and W v . Although tf u and * w 
have common eigenvectors, the ordering of the eigenvectors 
of ty v will in general be permuted relative to the ordering 
of the eigenvectors of van der Veen et al [4] have de- 
veloped a 2D extension of ESPRIT that exploits the fact 
that ¥ u and commute, since they have common eigen- 
vectors, to automatically pair the eigenvalues of * u with 
those of * v so that each resulting pair of eigenvalues is 
one-to-one associated with an incident wavefront thereby 
automatically pairing the direction cosine estimates with 
respect to the two array axes. 

Effectively, van der Veen et al propose a scheme for 
perturbing both 4 U and 4" so that the resulting per- 
turbed matrices commute. i~ l is then constructed from 
the eigenvectors of either perturbed matrix. Pre- and post- 
multiplying both the perturbed 4? u and the perturbed 4f v 
by *t and respectively, yields a pair of approximately 
diagonal matrices whose diagonal elements are automati- 
cally paired 1 . At this point, we will assume that the respec- 
tive eigenvalues of 4f u and 4t* have been correctly paired. 

The matrix whose columns are one-to-one related to the 
electromagnetic field vectors for each source denned by (10) 
may be estimated as 

Ai = i {E,t- 1 + E 2 i- 1 * u ' -hEat- 1 ^*} (15) 

Now, since the Jk-th column of *f ~* is the eigenvector asso- 
ciated with the k-th eigenvalue of 4f u , the k-th column 
of Ai above is associated with . Given the expression for 
the asymptotic form of A i in the far right hand side of (10), 
the Poynting vector for the fc-th source may be estimated 
via the vector cross-product between the top and bottom 
3x1 subvectors of the fc-th column of Ai according to (1). 
Normalizing the resulting cross-product vector to have unit 
length, the first and second components are high-variance 
but unambiguous estimates of the direction cosine of the 
fc-th source relative to the x and y axes, respectively. 

3.2. Disambiguation of ESPRIT Eigenvalues 

Let p = \p Sk ,Pv h ,Px k ] T denote the estimate of the normal- 
ized Poynting vector obtained by using the eigenvectors of 
the perturbed i M to decouple the signal eigenvectors com- 
prising Es and the subsequent column-wise vector cross 

1 In actuality, van der Veen et al employ a Schur decomposition 
as opposed to an EVD and there is much detail on how they 
approximately enforce commutativity that is not included here 
due to space limitations. 



product vector as described above. The first two compo- 
nents of this vector are used as coarse estimates to disam- 
biguate the candidate estimates of u* and vt described by 
(13) and (14), respectively. The disambiguated estimates 



a * = 7^** + roB T and *a = tAt^ + "°X 
2ttA A 2jtA A 

where m° and n° may be jointly estimated as 
Method 1: 



m ,n } = arg 



[• 



(16) 



where || || denotes the Frobenius norm, or estimated indi- 
vidually as 
Method 2: 

o def f . L A , . mAll 

m = |p t4 __^-_|j 

n " = [ min |p»* - 2^** - x|] (17) 

The search range for m and n in both cases is given by the 
far right hand side of (13) and (14), respectively. 

Method I tests up to a maximum of |2^ + lj 2 candidate 
pairs, whereas Method II only tests up to a maximum of 
2|2y + Ij candidates. Although Method II is computation- 
ally more efficient, it determines m° and n° in a separable 
manner and may thus yield "worst" estimates relative to 
the jointly based technique signified by Method I. However, 
simulations of both methods nave revealed no difference in 
their respective estimates when the inter- vector sensor spac- 
ing A is below a certain threshold. 

4. SIMULATIONS 

Simulation results presented in Figures 1 and 2 verify the 
efficacy of the extended aperture vector sensor array con- 
cept and attendant ESPRIT-based 2D angle estimation al- 
gorithm. The signal scenario involved two uncorrelated 
sources with the following location and polarization pa- 
rameters: $i = 41.4°, <f>i = 36.0°, 7i = 18° t m = -54°, 
$ 2 = 66.6°, fa = 46.8°, 72 = 72°,% = 126*. The SNR per 
component antenna was 20 dB for source 1 and 17 dB for 
source 2. For a given trial run, the 18 x 18 sample correla- 
tion matrix was formed from N = 32 snapshots. 

For each different value of the inter- vector sensor spacing, 
A, sample biases and sample variances for both u* and 0*, 
k — 1,2, were computed for each of the two sources from 
300 independent trial runs. For a given source, an RMS 
standard deviation was computed by taking the square 
root of the sum of the respective samples variances for u* 
and Oft. An RMS bias was also computed by taking the 
square root of the sum of the respective sample biases for 
Uk and Ok. Figure 1 shows the RMS standard deviation 
for source 1 as A is increased from a half- wavelength to 
over 50 wavelengths (100 half- wavelengths.) On a log-log 
scale it is observed that the RMS std. dev. decreases lin- 
early as the inter-vector sensor spacing is increased from a 
half-wavelength up to 32 wavelengths, at which point there 
is a factor of 80 reduction relative to the RMS std. dev. 
obtained with half- wavelength inter- vector sensor spacing. 
The RMS bias for source 1 plotted in Figure 2 exhibits a 
similar behavior and is an order of magnitude lower than 
the RMS std. dev. Note that the RMS std. dev. and RMS 
bias curves for source 2 are very nearly identical to that for 
source 1 but are not included here due to space limitations. 

Figures 1 and 2 also show the RMS standard deviation 
and RMS bias of the direction cosine estimates for source 
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1 obtained by the Poynting vector estimation procedure 
proposed in this paper. Both curves are observed to be 
relatively constant as the inter-vector sensor spacing is in- 
creased from A = A/2 to A = 32A. Both curves are 
observed to be well above the corresponding counterpart 
curves associated with the disambiguated direction cosine 
estimates, with the separation between the two widening as 
the inter- vector sensor spacing is increased from A = A/2 to 
A = 32A. This substantiates the claim that the Poynting 
vector estimation procedure yields high- variance but un- 
ambiguous direction cosine estimates in comparison to the 
direction cosine estimates extracted from TLS-ESPRIT's 
eigenvalues according to (13) and (14 J. Intuitively, the 
lower variance of the latter relative to the former is due to 
the fact that the direction cosine estimates extracted from 
the TLS-ESPRIT's eigenvalues depend on the size of the 
array aperture as measured by A: the larger A, the smaller 
the variance of //j k /2*rA. The relative high variance of 
the direction cosine estimates extracted from the Poynting 
vector estimates is intuitively due to the fact that they are 
inherently extracted from information provided by a single 
vector sensor which has no effective aperture. 

A breakdown phenomenon is observed at a thresh- 
old inter-vector sensor spacing of A — 32A (64 half- 
wavelengths.) The cause of this breakdown is still un- 
der investigation. The breakdown of the Poynting vector 
estimation procedure may be attributed to either fii = 

e J ^ Attl and pa = e J ^ Atta becoming nearly identical, or 

u\ = e**?*" 1 and i*i = e }3 P* v * becoming nearly identi- 
cal, when A is near 32 wavelengths. The breakdown of the 
disambiguation procedure at and above A = 32A may be 
attributed to picking the wrong grid value in a few of the 
300 trial runs. It is well known that the sample variance 
and sample bias estimators employed in these simulations 
are not robust and are highly sensitive to a few outliers. 

For purposes of comparison, simulations were conducted 
with an L-shaped array of 18 unpolarized scalar sensors, 
with 9 elements uniformly spaced at a half- wavelength along 
both the x and y axes. This results in an 18 X 1 array- 
manifold (same length as in the 3 vector sensor case) and 
an array aperture of roughly 8 half- wavelengths along ei- 
ther axis. Employing the same signal parameters presented 
previously, the unpolarized scalar sensor array yielded an 
RMS standard deviation of 5 x 10~ 4 for source 1. The 
3 vector sensor array and attendant algorithm achieved 
the same RMS standard deviation with A = 1. 5 A (3 half- 
wavelengths,) i.e., with less than half the aperture of the 
unpolarized L-shaped array! This performance gain is due 
to the polarization diversity inherent in a vector sensor ar- 
ray. The performance gain of the vector sensor array and 
attendant algorithm proposed in this paper relative to the 
L-shaped array with the same number of component anten- 
nas becomes even more dramatic as A is increased. Increas- 
ing the inter-element spacing in the latter case is not feasible 
since the ambiguities induced would not be resolvable. 

5. CONCLUSION 

A novel ESPRIT-based 2D angle estimation scheme was 
proposed involving the use of a right- triangular array of 
three vector sensors spaced much farther apart than a half- 
wavelength. Simulations were presented showing the excep- 
tional estimation accuracy one can achieve with a large ex- 
tended aperture allowed by the disambiguation facilitated 
by the electromagnetic field information provided by the 
vector sensors. Processing three vector sensors according to 
the method developed in this paper allows one to estimate 
the azimuth and elevation angle of up to a maximum of 
six sources. In order to handle more sources, the proposed 
scheme may be extended for a larger number of vector sen- 
sors spaced according to a rectangular array geometry with 



inter-vector spacings much greater than a half- wavelength. 
This will be presented in the journal version of the paper. 
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Figure 1: RMS standard deviation of tii and vi vs. inter- 
vector sensor spacing. 
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Abstract— We investigate direction-of-arrival (DOA) estimation 
involving nonuniform linear arrays, where the sensor positions 
may be noninteger values expressed in half-wavelength units, 
with some number of uncorrelated Gaussian sources that is 
greater than or equal to the number of sensors. We introduce 
an approach whereby the (noninteger) co-array is treated as the 
most appropriate virtual array when considering an augmented 
covariance matrix. Since such virtual arrays have an incomplete 
set of covariance lags, we discuss various completion philoso- 
phies to fill in the missing elements of the associated partially 
specified Hermitian covariance matrix. This process is followed 
by the application of an algorithm that searches for a specific 
number of plane wavefronts, yielding the minimum fitting error 
for the specified covariance lags in the neighborhood of the 
completion-initialized DOA estimates. In this way, we are able to 
resolve possible ambiguity and to achieve asymptotically optimal 
estimation accuracy. Numerical simulations of DOA estimation 
demonstrate a close proximity to the Cramer-Rao bound. 

Index Terms — Antenna arrays, array signal processing, direc- 
tion-of-arrival estimation, linear arrays, nonuniformly spaced ar- 
rays. 

I. Introduction 

CONSIDER an M-element noninteger-spaced nonuniform 
linear array, where the sensor positions may be noninteger 
values (measured in half-wavelength units). Since such an array 
generates up to |M(M — 1) distinct nonzero covariance lags, 
it has the potential [1] to estimate a superior number of uncor- 
related sources m, meaning that this number is greater than or 
equal to the number of antenna elements 

M<m< ±M(M-1). (1) 

The only currently existing techniques that may be applied to 
the superior case for such arrays are the multidimensional max- 
imum-likelihood search [2] and the geometry-invariant model- 
fitting idea [1], [3]. 

Here, we further investigate the geometry-based generalized 
augmentation approach (GAA) [4] and develop it for superior 
scenarios and noninteger linear arrays. We restrict our investi- 
gation to the case of uncorrelated Gaussian sources. Clearly, a 

Manuscript received January 15, 1997; revised October 6, 1999. This work 
was supported in part by a grant from INTAS SASPARC. Some results of this 
paper were presented at the ICASSP Conference, Munich, Germany, 1997. The 
associate editor coordinating the review of this paper and approving it for pub- 
lication was Editor-in-Chief Jos6 M. F. Moura. 

Y. I. Abramovich and N. K. Spencer are with the Cooperative Research Centre 
for Sensor Signal and Information Processing (CSSIP), Adelaide, Australia. 

A, Y. Gorokhov is with the Laboratoire des Signaux et Systemes, fecole 
Superieure d' Electric ite, Paris, France. 

Publisher Item Identifier S 1053-587X(00)02378-3. 



large class of practical problems is excluded by this assumption. 
However, some modeling restrictions are necessary to achieve 
superior-scenario capabilities. 

Another type of imposed restriction occurs for the existing 
high-order-moment algorithms [5], where the assumption of 
independency and non-Gaussianity is essential in achieving 
"blind" separation of such sources. It seems quite natural 
that although a conventional number of sources (defined as 
1 < m < M) can be localized, e.g., by MUSIC, under quite 
mild conditions, every attempt to extend localization capabili- 
ties beyond this limit involves some additional restrictions. 

The overall idea behind our approach is to use the noninteger 
co-array as the virtual array. Since the original set of specified 
(measurable) covariance lags is incomplete for this choice of 
virtual array, the DOA (spatial spectrum) estimation problem 
can be formulated in terms of a positive-definite (p.d.) Hermi- 
tian completion problem, analogously to the p.d. Toeplitz com- 
pletion problem for partially augmentable (integer) arrays [6]. 
Among possible completion criteria, we first investigate max- 
imum-entropy (ME) completion and explore the limitations of 
this approach. For the special class of noninteger arrays that are 
slight perturbations (in terms of antenna locations) of fully aug- 
mentable arrays, we introduce a better "minimum-deflective" 
completion. 

Since the proposed virtual array of co-array geometry is 
itself nonuniform, even with perfect completion, such an array 
can still suffer from manifold ambiguity [7]-[9]. Under our 
signal model, this means that for some particular scenarios, 
the co-array manifold (steering) vectors are linearly dependent. 
Therefore, MUSIC-type algorithms applied even to these 
perfectly completed covariance matrices will yield a set of 
ambiguous DOA estimates. 

Our goal in this paper is to find some technique of properly 
identifying the true sources for superior scenarios, provided that 
they are potentially identifiable. For a conventional number of 
(independent) sources, this question was addressed in [9]. Here, 
we extend that approach, whose main idea is to find the best 
single fit between the set of estimated covariance lags and the 
DOA and power estimates, using initial DOA estimates obtained 
by completion. In this way, we can determine those directions 
that actually correspond to the true signals, as opposed to those 
that have essentially zero power. Simultaneously, we are able to 
refine our DOA estimates to the asymptotically optimal bench- 
mark defined by the Cramer-Rao bound (CRB). 

We present results of numerical simulations that compare the 
actual accuracy achieved by our proposed algorithms against the 
CRB. For practical reasons, for any finite N and, hence, finite 
covariance lag estimation accuracy, we are also interested what 
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level of intersource correlation can be tolerated by our methods 
without a significant loss in DOA estimation accuracy. Corre- 
spondingly, the robustness of our basic assumption regarding 
source independency will be analyzed in a quantitative manner. 

This paper addresses the above-mentioned issues with the 
following structure. Background material is summarized in 
Section II. In Section III, we introduce our technique for ME 
p.d. completion. Section IV describes our "minimum-deflec- 
tive" completion as an alternative to ME completion. Section V 
is devoted to manifold ambiguity resolution and local DOA 
estimation refinement, as well as the nonidentifiability condi- 
tions that preclude such resolution by any means. Section VI 
presents the results of numerical simulations in comparison 
with the CRB, whereas Section VII summarizes the paper. 

II. Background 

Consider m narrowband planewave signals of power 
p = [pi , • • • , p m ] impinging on a nonuniform linear array 
of M identical omnidirectional sensors located at positions 
d = [0,^2, * • • ,c/m] measured in half-wavelength units. We 
assume that the number of plane-wave sources m is known. 
Adopting the commonly used data model [10], we have 

y(t) = S($)x{t)+ V (t), for < = 1, • • • ,iV (2) 

where 

x{t) = MO, • • • ,*m(0] r , v(t) = ■ • • :2/m(0] t 

v(t)=[vi(t)r-',VM(t)] T (3) 

x j(t) U = 1, * * • , m ) is the complex signal amplitude of the 
jth plane wave, and where yu(t) and rfk(t) (k = 1, • • • , M) are 
the sensor output and the noise at the kth sensor, respectively. 
To permit DOA estimation in the superior case (m > M), we 
restrict ourselves to the class of independent Gaussian signal 
amplitudes x(t) e C mXl such that 

We assume that the additive noise rj(t) € C Mxl is white and 
Gaussian 

eWMfm-ff"- ( 5) 

The array manifold matrix is S(0) = [s($i ),♦••, s(0 m )] € 
C Mxm , where each constituent "steering vector" s(0k) is de- 
fined as 

= I 1 ' expfMHfctiy), • • • , exp(z7rd M ™i)] T (6) 

with w - sin0 € [-1,11- 

According to this model, the M-vafiate spatial covariance 
matrix 

R = SPS H +p 0 I M (7) 

is p.d. Hermitian. Note that in our case of m > M, the noise-free 
covariance matrix SPS H is generally of full rank. Given N 
independent samples ("snapshots"), the sufficient statistic for 
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DOA estimation is the M-variate direct data covariance (DDC) 
matrix 

*=Ji E »(*)»"(*)• ( g ) 

Under our assumptions, R is characterized by a complex 
Wishart distribution NR ~ CW(iV, Af , R), nondegenerate for 
N > M, and thus, the CRB canbe readily calculated [10] to 
provide the asymptotic bound ori^^^^pS^^^^is^Yacy for 
every "locally identifiable" scenario (i.e. with nondegenerate 
Fisher information matrix) [9]. We use the CRB as our DOA 
estimation accuracy benchmark. 

To illustrate the problems that arise in the superior case, con- 
sider the four-element noninteger array introduced by Proukakis 
and Manikas [7] 

d 1 = [0,1.2,3.4,4.6]. (9) 

In considering the superior scenario, we begin with the obser- 
vation that the co-array of di is Ci = [0, 1.2, 2.2, 3.4, 4.6], i.e., 
this geometry specifies the set of L = 5 covariance lags 

*" = [ r j]jev = [r( Ci )]j=i, -,l 

— I r (0) 5 r (1.2) j r (2.2) 5 r (3.4) 5 r (4.6)] ( 1 0) 

'(where L is the number of co-array elements and, of course, is 
equivalent here to M a ), which is the number of virtual array 
sensors, which means that the maximum number of potentially 
identifiable sources in this case is m max = 4, whereas the 
number of (real) degrees of freedom is 2L — 1 = 9 [1]. In addi- 
tion, note that the co-array geometry ci still meets the sufficient 
condition for manifold ambiguity defined in [7] 

c L > L-l. (11) 

The Proukakis-Manikas technique for determining ambiguity 
generator sets (AGS's) may be applied to the co-array to find 
that there is a single AGS associated with ci, which is 

wags = [-0.8391, -0.4043, 0.0304, 0.4652, 0.9]. (12) 

Thus, for the (virtual) array cu any four DOA's from the 
set wags result in all five MUSIC peaks, even if the true 
five-variate covariance matrix is used. While the array di 
and its co-array ci share the same AGS wags, there is an 
important distinction; the ambiguity rank is equal to three for 
the array and four for the co-array. This means that any three 
manifoldly ambiguous independent sources obtained via the 
Proukakis-Manikas technique are identifiable fordi. 

For the superior case of four sources with d u we have to first 
devise a new DOA estimation approach, given the set of covari- 
ance lag estimates f, and, second, to explore avenues of ambi- 
guity resolution given a manifoldly ambiguous scenario such as 

ti>4 = [-0.8391, -0.4043, 0.0304, 0.4652]. (13) 

These two problems can be tackled independently since if the 
full set of ambiguous DOA's (five in this case) is somehow iden- 
tified, all that remains is to decide which of them are true and 



ABRAMOVICH et ai: DOA ESTIMATION FOR NONINTEGER LINEAR ANTENNA ARRAYS 



945 



which are spurious. For an identifiable scenario, the manifold 
ambiguity resolution problem will have a unique solution. 

In the next section, we discuss an approach capable of (pos- 
sibly ambiguous) DOA estimation for a superior number of un- 
correlated sources. We deal with manifold ambiguity resolution 
in Section V. 

III. Maximum-Entropy Positive-Definite Hermitian 
Completion 

When the co-array c is selected as the virtual array, the set 
of specified covariance lags r defined by the array d is incom- 
plete for the virtual array, as is the L-variate augmented Hermi- 
tian covariance matrix H . For our example noninteger geometry 
di = [0, 1.2,3.4,4.6], the partially specified five-variate aug- 
mented covariance matrix is 



H = 



r (o) 

r (l.2) 
r p.2) 

r (4.6) 



r U-2) 
r (0) 
? 

r {2.2) 
r (3.4) 



r (2.2) 
? 

'(1.2) 



T '(3.4) 
r (2-2) 
r (1.2) 

r (1.2) 



r (4.6) 

r (3-4) 
? 

r (12) 
r (0) 



(14) 



since the covariance lags H 2) 3 = r( 1>0 ) and #3,5 = r (2 .4) 
are not presented by the array d\. Thus, DOA estimation for 
noninteger NLA's is closely connected with the mathematical 
completion problem for partially specified Hermitian matrices. 
Since maximum likelihood (ML) completion is not directly at- 
tainable, we will instead propose some other ad hoc completion, 
followed by ML refinement. 

We begin with maximum-entropy (ME) completion, 
although, in general, as a result of ME completion, we may ob- 
tain a solution that cannot possibly correspond to a plane-wave 
model. Moreover, since entropy is defined by the determinant of 
the covariance matrix for a Gaussian process, ME completion 
may significantly increase the effective signal subspace, i.e., 
the number of eigenvalues significantly greater than the white 
noise power po or, on the contrary, may increase the virtual 
"noise floor" in the completed matrix. This is why we have 
to be cautious applying ME completion, and we should be 
prepared to treat it only as an initial step. 

Let the general virtual array cf be specified by the virtual 
sensor positions (j = 1, • * • , M'); then, the set of all possible 
p.d. Hermitian completions H may be written as 



H = \z: H(z) = H 0 + £ (Re H pq E™ 
P,qes 

P<9 



+ i Im H M E P J ) > 0 



(15) 



where 



z — 



Re Hr, 



lmH„ 



PV J pq€Sp<q 



and where e p = [0, • • • , 0, 1, 0, • ■ • , 0] is the M'-variate basis 
vector with a unit entry in the pth position, Ho is the augmented 
covariance matrix with zeros in place of the unspecified ele- 
ments of H, $ is the set of specified elements {p, q) of if, and«S 
is the set of unspecified elements in the initial incomplete aug- 
mented covariance matrix H. The specified elements of H are 
defined by the standard direct augmentation approach (DAA) 
[11] operating on the M-variate covariance matrix R (7) 
M 

M pqes - — • 



(17) 



Y, tivp-Vqidj-dk) 



j,k=l 



In practical situations, we apply the DAA to the DDC matrix R 
(8) to obtain H. 

Suppose we label each of the missing lags (j>q € <S;p < 
q) from 1 to £, which is the total number of missing lags. In 
our example, we then have the set of £ = 2 missing lags S — 
{#2,3 = r (i.o) > #3,5 = r (2.4)}* Then, we may introduce the 
following notation instead of (15): 



,.( 



21 



= ^ z: H{z) = H 0 + Ys Z 3 F > > 0 



(18) 



where 



*3 ~ \ iE P1 



Re V pqeS;p<q> 

Im r pqes,p<q' 



(19) 



for j =£+l,--',2£ 

for j = 1,...,/ 
fori = / + l,---,2^. 

Let us assume for the moment that the feasibility condition for 
the p.d. completion of H(z) is satisfied, i.e., there exists at 
least one solution z for the linear matrix inequality (LMI) [12] 
H(z) > 0. For exact covariance lags (deterministic comple- 
tion), this is always true since there exists the p.d. completion 
corresponding to the true DOA's and the virtual array. When the 
feasibility condition is guaranteed, the following theorem may 
be immediately applied. 

Theorem 1: [13], [14] Let if be a partially specified posi- 
tive-definite matrix, all of whose diagonal elements are spec- 
ified, and suppose that H admits at least one positive-definite 
completion. Then, among all such completions, there is a unique 
one with maximum determinant. Furthermore, this determinant- 
maximizing positive-definite completion Hm e is characterized 
by having a zero in its inverse in every position where H has an 
unspecified element. 

This theorem generalizes the well-known property of the 
Burg ME spectrum, whereby an ME-extended Toeplitz co- 
variance matrix has zeroes in its inverse for every unspecified 
subdiagonal (the auto-regressive property) [15]. 
According to Theorem 1 , the function 

det H~ l (z) } for z eH 
otherwise 



<P(z) 



= (log 
\oo, 



(20) 



is strictly convex on the feasible set and therefore has the unique 
minimizer 



(16) 



z M e = arg min <p(z) = axg max det H(z) 



(21) 
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which is called the analytic center of the LMI H (z) > 0. 
For Gaussian distributions, we may evidently treat the analytic 
center as the ME completion since the entropy for Gaussian-dis- 
tributed random variables of zero mean is given by [15] E = 
| log(detif), where H is the correlation matrix. 

Let us now turn to the problem of computing the analytic 
center of the LMI. Since the feasibility conditions are assumed 
to hold at this stage, the search for the analytical center can be 
split into two subproblems. First, define an arbitrary feasible 
point z such that H(z) > 0. Second, given z, any convergent 
procedure that minimizes <f>(z) (taking into account the feasi- 
bility condition) will converge to the analytical center zme- For 
details of the Ellipsoid and Newton algorithms used to solve 
these two convex programming problems, see [12] and [16], 
but to summarize, convex programming solves the generalized 
problem of finding the ME completion for an arbitrary feasible 
incomplete covariance matrix. 

When the sufficient statistic is provided in the form of a 
M-variate sample covariance matrix R characterized by a 
complex Wishart distribution CW(N, M, R), the set of sample 
specified covariance lags f is usually obtained via the direct 
augmentation approach of (17) with sample covariance lags 
instead of exact values. In this case, the feasibility condition 
for the existence of a p.d. Hermitian completion is no longer 
guaranteed to hold. 

The question of which partially specified Hermitian matrices 
may be p.d. completed is addressed in [14]. It was shown that 
if the diagonal elements are specified, if the principal minors of 
the specified elements are positive, and if the undirected graph 
of the specified elements is chordal, then a p.d. completion ex- 
ists. This last condition always seems to be satisfied for the aug- 
mentation approach simply because the p.d. completion always 
exists for the exact covariance lags (with the same undirected 
graph, of course). Instead of forming all principal minors, we 
may simply associate the success of the Ellipsoid Algorithm in 
finding a feasible point z with the feasibility condition men- 
tioned above. 

Hence, if the necessary and sufficient condition fails (the fea- 
sible point z is not achievable by the Ellipsoid Algorithm), we 
need to modify the set of sample specified covariance lags in 
order to achieve feasibility with the minimum possible deflec- 
tion from the initial set of (maximum likelihood) estimates f. 
Several approaches for such a modification may be proposed. 
For example, by introducing 

H(z) = Ho + Yl z '& + E *i p i < 22 > 
jes j( =s 

similarly to (18), we may find the minimum-deflective point by 
the solution of the optimization problem 

find min V) (zj) 2 subject to H(z) > 0. (23) 

The Ellipsoid Algorithm provides a straightforward approach to 
find the unique optimum solution. 

Another possibility is to define the minimal diagonal loading 
for the matrix H(z) to achieve feasibility, which is also a 
standard convex optimization problem. In general, all these 
approaches may provide different solutions in particular cases. 



Nevertheless, we do not expect statistically different results 
in the asymptotic domain, where in accordance with the 
Cramer-Rao bounds, we expect reasonable convergence to the 
true DOA's. 

We note that this ME Hermitian completion approach is valid 
for arbitrary-geometry arrays, e.g., 2-D. 

IV. ME and MD Positive-Definite Hermitian Completion 

The noninteger NLA di = [0, 1.2,3.4,4.6] introduced by 
Proukakis and Manikas [7] is an example of a poor geometry 
choice for DOA estimation purposes since it has only five 
distinct covariance lags rather than the maximum possible of 
\M{M - 1) + 1 = 7. In order to increase the number of 
potentially identifiable sources, all redundancies should be 
avoided. 

Let us instead turn our attention to the noninteger array 

d 2 = [0,1.09,3.96,5.93] (24) 

which is recognizable as a slightly perturbed version of the "op- 
timal-lag" (nonredundant and fully augmentable) integer array 
d = [0, 1,4,6] investigated by Moffet [17]. In practical appli- 
cations, such geometry perturbations are usually caused by an- 
tenna positioning errors. This quasi-fully augmentable array d 2 
is also nonredundant and generates the full set of seven lags: 

c 2 = [0,1.09,1.97,2.87,3.96,4.84,5.93]. (25) 

Naturally, this co-array is a slightly perturbed version of 
the seven-element ULA. In this case, the array manifold 
ambiguity condition d,M > M — 1 is satisfied, whereas the 
co-array manifold ambiguity condition cl > L — 1 is not 
satisfied, and therefore, the only currently existing method (the 
Proukakis-Manikas technique) is unable to find any AGS's for 
the co-array c 2 . 

The class of noninteger linear arrays that do not satisfy the 
co-array manifold ambiguity condition may be called GIMP 
identifiable since, like fully augmentable (integer) arrays, their 
co-arrays appear to have no AGS's, and therefore, they would 
be globally identifiable (identifiable regardless of the DOA set). 
We have conditioned this assertion by the acronym GIMP for 
two reasons. First, because this is true assuming Gaussian inde- 
pendent sources, and second, because the Proukakis-Manikas 
technique does not yet guarantee to produce an exhaustive list 
of AGS's, i.e., the manifold ambiguity condition d M > M - 1 
has been proven to be sufficient, but it is not known whether it 
is necessary. 

While the potential capabilities of the array d 2 are obviously 
more attractive than those of rfi, the problem of L-variate aug- 
mented covariance matrix completion is slightly different in this 
case. Nonredundant arrays have the number of co-array ele- 
ments L = ~M(M — 1) + 1, whereas the number of unspecified 
covariance lags in the L-variate augmented covariance matrix is 
t = ^(L—1)(L — 2). Even prior to computational examples, we 
may expect for M >- 1 that the extreme incompleteness of the 
augmented matrix might lead to an ME completion that is quite 
inconsistent with the DOA estimation problem. In fact, the in- 
verse of the ME completion (H]^ E ) will have exactly £ zeroes 
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in the position of each missing lag and is obviously radically 
different from the structure of the inverse of the true covariance 
matrix H~xact f° r m independent plane waves. 

In view of these misgivings regarding the appropriateness of 
the ME completion, it seems quite natural to exploit the close 
proximity of the perturbed array geometry to the corresponding 
integer array. In the case of quasifully augmentable arrays, the 
following minimum-deflective (MD) completion treats the true 
Hermitian matrix H exac t as a slightly perturbed Toeplitz matrix. 
We can write each missing lag of the augmented Hermitian ma- 
trix as a perturbation such as 

H pq = H jk + z n + iz n + £ for n,p,q€S; j,keS (26) 

where the nearest specified covariance lag is used 



mm 

jk 



\cj - c*| - \c p - c q \ 



(27) 



In fact, the same approximation may be made for any type of 
noninteger array, such as di, where 

#2,3 = r(i.Q) = ?'(1.2) + zi + ™2 (28) 

#3,5 = **(2.4) = r (2 . 2 ) + z 3 + iz 4 . (29) 

In general, within the set of specified covariance lags «S, we can 
always find the "nearest" existing lag to any given missing lag 
of the virtual array (the co-array) and use this completion H T 
instead of H$. Naturally, this completion may not be a feasible 
covariance matrix since this process cannot guarantee the posi- 
tive-definiteness of Ht- 

The definition of the minimum-deflective solution zmd that 
makes the L-variate Hermitian matrix H(z M d) = Hmd non- 
negative definite is 

find min ^ |jgj| 2 subject to 



2i 



(30) 



H(z) = H T + Y, ZjFj > 0- 

For our example array, d 2 ,£ = \{L — 1)(L — 2) = 15. 

It is now clear that we may treat this MD completion as a 
Pisarenko-type completion [6]. Formally, the problem (30) is 
identical to the minimum-deflective initial point search with the 
initial matrix H 0 replaced by H T . We may make such a com- 
pletion more effective by modifying the diagonal element by a 
fixed correction in noise power (ro — Po), or by an optimised 
deflection (f 0 — A), provided that the feasibility condition is not 
violated. 

To summarize, convex programming techniques provide 
computationally effective completions of partially specified 
Hermitian covariance matrices of the so-called virtual array 
that may be used to initialise DOA estimation. 

V. Manifold Ambiguity Resolution and Local DOA 
Refinement 

There are four major reasons why DOA estimates obtained 
by the above completion methods need to be refined. First, we 
have shown that the virtual array can suffer from manifold am- 
biguity, even for the true (deterministic) L-variate covariance 



matrix. Second, even for the true covariance lags r specified 
by d, the ME- and MD-completion techniques do not neces- 
sarily yield the true DOA estimates, i.e., these techniques lead 
to biased estimates. Third, ME completion increases the deter- 
minant of the completed matrix H ; in turn, this has the potential 
to introduce spurious sources. This "self-induced" ambiguity is 
a more serious problem than co-array manifold ambiguity since 
the latter exists only for very specific DOA sets (defined by the 
AGS's), whereas the former phenomenon can arise in an arbi- 
trary DOA scenario. 

The final reason is that the DAA adopted for estimating spec- 
ified covariance lags has been proven to be far from asymptot- 
ically optimal [11], even in its pure form. Naturally, the result 
of the additional step of completing a partially specified ma- 
trix cannot generally improve the overall DAA performance. 
This qualification is important since the biased nature of DOA 
estimates obtained via completion may, in principle, lead to 
super-efficient estimates [18] in particular cases. 

Given the set of initial DOA estimates obtained by a com- 
pletion method, we next have to solve two problems. First, 
we need to be able to resolve either (or both) of the ambi- 
guity types mentioned above (co-array manifold ambiguity 
and self-induced ambiguity). Second, having eliminated the 
spurious DOA estimates, we need to increase the accuracy of 
the properly identified estimates in order to reduce any bias in- 
troduced by completion and to reduce the stochastic estimation 
errors as close as possible to the CRB limit. By assuming in 
this paper that the CRB is finite, we are deliberately excluding 
situations where the Fisher information matrix is degenerate, 
i.e., where the scenario is locally nonidentifiable [9]. 

A. Manifold Ambiguity versus Nonidentifiability 

Here, we consider the problem of determining whether any 
manifoldly ambiguous scenario (with independent Gaussian 
sources) is, in fact, nonidentifiable. We will shortly see that 
the approach presented here is also capable of resolving the 
"self-induced" ambiguity created by ME completion. 

According to the model described by (4) and (7), the set of 
specified spatial covariance lags may be presented in the form 

r f = A(0)p forp > 0 (31) 

where r' = [r 0 - po> r\ > * * * > ^l-i] t is the white-noise modi- 
fied L-variate vector of specified covariance lags, and A(0) = 
[a(#i), • • • , a(0 m )] is the L x m co-array manifold matrix, with 

a(0j) = [1, exp(i7rc 2 Wj), • • • , exp(i7TCLWj)] T . (32) 

Pointwise nonidentifiability therefore occurs at the point 0+ if 
there exists a different set of m DOA's 0**^0+ and two sets of 
(possibly different) non-negative powers p+ such that 

A{9*)p. = A(0 f M. (33) 

Obviously, this will occur when there exist the pairs 0+ ^ 0+ and 
such that their corresponding multivariate distributions 
are identical, i.e., 

f(y\o.,p.) = MK,A) 04) 



since 



(35) 
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Hence, (33), which is identical to (35) here, is the necessary and 
sufficient condition for nonidentifiability. 

Note that we are discussing the case where the pair 0* , 0* may 
be chosen from a set of measure zero, where the Kullback di- 
rected divergence (KDD) 7(0*, 0*) between the corresponding 
Gaussian probability densities is strictly equal to zero [ 1 9] , [20] , 
and where the "generalized ambiguity" [20] is strictly equal to 
one. 

Let us illustrate the distinction between the co-array manifold 
ambiguity condition (11) and the nonidentifiability condition 
(33) starting from the simple case when the set 0^ differs from 
0* by a single DOA so that the set 0* = 0* U 0* consists 
of (ra + 1) DOA's. In this case, the necessary and sufficient 
condition for nonidentifiability (33) may be rewritten as 



A{h)h = 0 



where 



h = 



P+i 
p*2 - At 



P*m " P*m-l 
-P*m 



(37) 



with J?*,!** > 0. This condition can be rephrased as occurring 
when the complex homogeneous linear system (36) admits a 
solution h that is real and has the property 

hhm+i < 0. (38) 

It is now clear that our application of the manifold ambiguity 
condition [7] to the co-array 

rank A(0*) < m (39) 

that ensures the existence of an arbitrary nontrivial complex so- 
lution to (36) is only a necessary condition for nonidentifiability 
since it does not guarantee the extra two properties of h. Thus, 
the Proukakis-Manikas technique applied to the co-array will 
generate only potentially nonidentifiable scenarios. 

Since the solution h is required to be real, (36) can be 
rewritten as 



A(0+)h = 0 



where 



A{K) = 



Re A(0*) 
Im A{K) 



e ^(2L-l)x(m+l) 



(40) 



(41) 



(the first row of the imaginary part is dropped since it is trivially 
equal to zero). An arbitrary nontrivial real solution of the new 
system (40) exists if and only if 

rank .4(0*) < m. (42) 

This condition is clearly more stringent than (39); in fact, all 
AGS's generated for the co-array geometry examples men- 
tioned below did not meet condition (42), whereas obviously, 
they did satisfy (39). Hence, (42) can be considered to be an 
additional necessary condition for nonidentifiability, and the 
list of co-array manifoldly ambiguous scenarios generated by 
the Proukakis-Manikas method may now be checked against 
(42) to further select potentially nonidentifiable scenarios. 



Suppose that ,4(0*) is tt-rank-deficient, i.e., rank .4(0*) = 
(m + 1) - k, with k > 1, so that every possible (real) solution 
to (40) can be written as 



h = $x 



(43) 



where $ 6 7£( m+1 ) x * is the fundamental solution to the system 
(40), and x € K KXl . The purpose of this step is to allow us 
to search for a possible solution h of (40) with the two extra 
constraints by formulating a linear programming (LP) problem. 
Using a standard "trick," we present the arbitrary real vector x 
in the form x = xi - X2 for some x\,x 2 € 7£+ xl (i.e. posi- 
tive-valued vectors). We can now verify that the following LP 
problem is appropriate for solving $x - h = 0 with /ii/i m +i < 
0 



(36) Find min (A x + A 2 ) for - x 2 ) 



0 




0 




2/3 














0 




0 


0 




0 




-2/4 




A 2 



= 0 (44) 



where 2/1,2/2 € T^ 7 * -2 ^* 1 and 2/3,^4, Ai,A 2 > 0 since if and 
only if the minimum solution is precisely zero will the solution 
h be of the required structure, namely 



h = $x = $(xi - X2) = 



2/3 
V2 -Vi 
-2/4 



(45) 



where all the y's are positive valued. Note that this LP problem 
is guaranteed to have a solution since for nonzero Ai , A 2 , we can 
always write <J>(xi - x 2 ) = zi - z 2 for some z x ,z 2 € 7£ mxl . 

To summarize, the existence of a zero solution to the LP 
problem (44) is the necessary and sufficient condition of non- 
identifiability for the co-array manifoldly ambiguous scenario 
0* =0* U 0*. 

Note that in this same situation but with the different a priori 
model where the powers p+j (j = 1, * • • , m) must be identical 
(p+j = p say), nonidentifiability occurs under the more stringent 
condition 



fci+Am+i=0; hj=0 fori = 2,. 



, m 



(46) 



since the solution could be normalized such that hi = p > 
0. The significant difference between the two conditions (38) 
and (46) is due to the fact already discussed in [9], [21], [22]: 
that identifiability conditions strongly depend on the statistical 
signal model and are not defined by the array geometry alone. 

Up to this point, we assumed that 0* = 0* U 0^ consisted of 
(m + 1) DOA's. In general, it will consist of m + 1 < £ < 2m 
DOA's so that only (2m - £) DOA's belong to 0* n 0* . For the 
arbitrary power model, the rank-deficiency of .4(0*), k may be 
greater than one: 



rank-4(0*) = £-«; 



(47) 



for 1 < k < t — 2, whereas the general necessary and sufficient 
condition for nonidentifiability is 

0 



$x 



+ 



= 0 



(48) 



where $ 6 H tyiK is the fundamental solution to the system (40), 
x € ll KX \ andp +3 pi e 7l+ xi . As before, the solution to (48) 
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may be found by solving the following linear programming (LP) 
problem: 



t-m 



Find min ^ (Ay + A 2j ) for <&{xi - x 2 ) 

3=1 





" 0 " 




" o " 








' A x ' 


+ 


yi 






+ 


0 


+ 


0 




0 




0 








-A 2 



= 0 (49) 



where xi,X2 € 



Ai, A 2 G ^ m)xl . The nonidentifiability condition (40) is 
satisfied if this minimum is zero. Thus, for any locally identi- 
fiable scenario 0+ with arbitrary powers p+ and rank-deficient 
matrix A(0+), pointwise nonidentifiability occurs if and only 
if the solution of this LP problem is equal to zero. In other 
words, the necessary and sufficient condition for a locally 
identifiable uncorrelated Gaussian superior scenario to be 
pointwise nonidentifiable is 



Y, ( A i; + A 2i) = o. 



(50) 



Any such scenario that does not meet this condition is pointwise 
identifiable. 

Note that if we had a method for exhaustively listing all 
AGS's that satisfy (39) for a given co-array, then deciding 
whether an array is always identifiable or not would reduce 
to testing each AGS against the condition (42) and then the 
LP zero solution. We have already mentioned that the only 
existing technique for calculating AGS's [7] does not guarantee 
finding the full set of all possible AGS's. Since our study relies 
on this Proukakis-Manikas technique, complete identifiability 
of any given array can only be given in terms of the same 
(possible) subset of AGS's. However, any given array and 
DOA scenario can be easily tested. For example, the co-array 
geometry c\ results in the single AGS wags — sin0* by 
the Proukakis-Manikas technique, but this fails even to meet 
the second necessary condition (42). This means that the 
four-source scenario w± is globally identifiable, despite its 
co-array manifold ambiguity. 

B. Co-Array Manifold Ambiguity Resolution and Local DOA 
Refinement 

Given any co-array manifoldly ambiguous (but identifiable) 
scenario, we now wish to find an algorithm that can identify 
the true m DOA's 0+ amongst the set of £ > m + 1 DOA's 
0+ provided (for example) by MUSIC applied to the properly 
completed augmented covariance matrix H, 

Since, for identifiable scenarios, there is a unique solution to 
the equation r' = A(0+)p+, a straightforward idea is to consider 
the equation r' = A(0+)p+ and note that exactly (£ — m) com- 
ponents of the ^-variate positive- valued vector p+ will be equal 
to zero. Naturally, in the practical case when r and 0* are es- 
timated quantities, we should search for the best fit over all £ 
elements ofp % and associate the m greatest values with the true 
source powers. 



Instead of the least-squares optimization introduced in [23], 
we suggest that this "diagonal fit" problem (DFP) be solved 
directly by another LP 



where 



Find min/ x subject to Alpx = r f r and 



(51) 



0 T 



0 T 
-h(L-i)\ 



Ref' 
Imf 



6 ft(M-«xl 



/ = 



I 4L-4 



e 7£(2L-l)x(*+4Z,-4) 

(52) 
(53) 

(54) 



Both Im A(0*) and Im r consist of (L - 1) rows since we have 
deleted the row trivially equal to zero. Since the first £ elements 
of x correspond to the powers, as the number of snapshots N 
(composing the DDC matrix R) tends to infinity, precisely k of 
these elements will tend to zero, whereas the remaining m of 
these elements approach the powers corresponding to the true 
DOA's. Thus, the solution of the DFP is obtained by choosing 
the m greatest values amongst the first £ elements of x. 

This approach was inspired, in part, by Fuchs' model-fitting 
method [3], [24], as well as some of our early investigations 
into LP applications for array processing [25], This basic LP 
problem (51)-( 54 ) can be modified in different ways; for ex- 
ample, in addition to the mismatch minimization, we could also 
minimize the smallest element within the first £ elements of x. 
It should be clear that the solution of the DFP does not depend 
on the array geometry and could be easily generalized for arbi- 
trary nonlinear arrays (circular, 2-D, etc.), provided that identi- 
fiability is guaranteed for the given scenario. 

For the superior case, where the DOA estimates are initial- 
ized by a completion technique, ambiguity resolution should 
be followed by some sort of ML refinement to achieve near 
asymptotic efficiency for the final DOA estimates. The algo- 
rithm shown at the bottom of the next page includes such a local 
refinement (at Step 5). 

VI. Simulation Results 

To demonstrate the efficiency of the proposed techniques, 
consider the manifoldly ambiguous (but identifiable) scenario 
[7] 

(*! =[0,1.2,3.4,4.6] 

w 4 = [-0.8391, -0.4043, 0.0304, 0.4652]. (55) 

The problem is to estimate the true four sources tu 4 
(m = M) amongst the set of five ambiguous directions 
given by the co-array AGS wags = [-0.8391,-0.4043, 
0.0304,0.4652,0.9]. 

First, we investigate deterministic ME completion for this 
scenario at a SNR of 20 dB common to all four sources. The 
dotted line in Fig. 1 shows the ME spectrum of the augmented 
covariance matrix completed by the maximum-entropy criterion 
H M e- For comparison, the dashed line shows M£(i/ e xact), 
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ME, d = [0.0, 1 .2, 3.4, 4.6], 20dB 
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Fig. 1. ME spectra of the exact co-array covariance matrix H exact and the 
ME-completed matrix Hme for the AGS associated with the co-array of d. 
Arrows indicate the four true DOA's. 



which is the ME spectrum of the true covariance matrix, calcu- 
lated using the co-array (our virtual array) ci . In this case, the 
ME-completed and the true ME spectra are essentially identical. 
Both spectra exhibit five strong peaks in the expected directions 
wags- The corresponding MUSIC pseudo-spectra (which is not 
shown) for H exact and Hme (and, indeed, the minimum-deflec- 
tive completion Hmd) are all practically identical. The simple 
conclusion is that for this scenario, ME completion of the two 
missing lags produces results consistent with the DOA estima- 
tion problem. 

To emphasize the point that ME completion is an ad hoc cri- 
terion, not necessarily consistent, and prone to "self-induced" 
ambiguity, we next consider the same scenario but with three 
DOA's: 



di =[0,1.2,3.4,4.6] 

ti/ 3 = [-0.8391, -0.4043, 0.0304]. 




(56) 



-1 -0.75-0.5-025 0 0.25 0.5 0.75 1 
DOA (w) 

Fig. 2. Same as Fig. 1, except involving the first three DOA's in the AGS 
u>ags (12). 



Recall that the ambiguity rank of the AGS wags is equal to 
three for di and four for ci, i.e., this scenario is manifoldly 
ambiguous for the array but not for the co-array. Nevertheless, 
ME(Hme) (as shown in dotted line in Fig. 2) once again 
exhibits all five peaks, unlike the ME spectrum of the true 
covariance matrix ME(H exact ) (dashed line). Meanwhile, the 
MUSIC pseudo-spectrum of Hme assuming m = 3 sources 
shows four strong peaks (which is not shown), unlike the 
three true DOA's indicated by MUSIC (H exact ). Moreover, 
MUSIC (Hme), assuming four sources, presents all five peaks 
in the directions wags, similarly to MUSIC (Hme) in the 
previous four-source scenario. 

Thus, maximization of the determinant of the partially speci- 
fied covariance matrix has the effect of changing this potentially 
co-array-unambiguous three-source scenario into an ambiguous 
four-source scenario, with the true directions being a subset of 
some new ambiguous DOA set. 



Ambiguity Resolution with Local Refinement (ARLR) 
Step 1: Given the ME- or MD-completed DAA-augmented covariance matrix H(17) obtained from the 
DDC matrix R(8), calculate the initial set of£>m DOA estimates 0* by applying MUSIC to the 
L -variate matrix H and selecting all local maxima from the pseudo-spectrum. (If £ = m, then there is no 
ambiguity: stop.) 

Step 2: Solve the LP problem (51)— (54) to find the best-fit source power estimates p*. 

Step 3: Select the m greatest elements of p + and treat the corresponding m elements of 0* as the initial 

set of unambiguous DOA estimates 

Step 4: Calculate the set of refined unambiguous source powers p^ 0 ^ by solving the LP problem (51)-(54) 

with the m unambiguous DOA's 0 (o) substituted for the £ ambiguous DOA's 0*. 

Step 5: Apply the local ML search in the neighborhood of ($^°\p^) to obtain as described for 

example in [1 1] (therein referred to as Step 4 of the ATML algorithm). 
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TABLE I 

DOA AND POWER ESTIMATES OBTAINED 

via ME and MD Completion for Exact Covarjance Lags 
for the Fig. 1 Scenario 



MUSIC, d= [0.0, 1.2. 3.4, 4.6], 20dB 





& MD 


a ME 












W m 


P. 


p* 


P* 


-0.8391 


-0.8348 


-0.8392 


100.0000 


74.2378 


100.0004 


-0.4043 


-0.4096 


-0.4044 


100.0000 


96.6244 


100.0101 


0.0304 


0.0356 


0.0304 


100.0000 


100.0437 


99.9921 


0.4652 


0.4608 


0.4652 


100.0000 


103.1554 


99.9974 


0.9000 


0.9044 


0.9000 


0.0000 


25.9388 


l.OxlO" 10 



Interestingly, in this four-source experiment, both ME and 
MD completions retain the true eigenspectrum quite accurately. 
In contrast, ME completion of H in the three-source scenario 
results in a very significantly increased fourth eigenvalue (as a 
consequence of determinant-maximizing) and has allowed ap- 
preciable energy to leak into the noise subspace, thus creating a 
"self- induced" spurious extra source. Even MD completion dis- 
plays a similar phenomenon, although not to the same extent. 
This proves our assertion that ME (and MD) completion is not 
always consistent with DOA estimation, and both should be re- 
garded as ad hoc criteria. 

Note that the three-source scenario is identifiable with re- 
spect to di ; thus, all spurious directions can be easily removed 
by the "diagonal-fitting" method; more details on manifold am- 
biguity resolution in the conventional case can be found in [9]. 

Now, let us consider the accuracy of DOA estimation by our 
techniques. Table I presents the DOA and power estimates ob- 
tained by ME and MD completion using the true (determin- 
istic) covariance lags and the diagonal-fit method for the above 
four-source scenario. In this particular case, ME completion 
yields almost ideal results, whereas MD completion is less ac- 
curate. However, the worst DOA error under the MD scheme 
(0.0052) is comparatively small since it corresponds to the CRB 
for almost 500 snapshots. Note that the inferior MD-compIeted 
DOA estimates have resulted in a considerable redistribution of 
power from the first (true) DOA to the fifth (false) DOA. 

In order to demonstrate that the relative merits of ME and 
MD completion are scenario-dependent, we next consider the 
arbitrarily chosen unambiguous four-source scenario 

di = [0,1.2,3.4,4.6], t» 4a = [-0.8,-0.3,0.1,0.4]. (57) 

Fig. 3 illustrates the MUSIC pseudo-spectrum of H exact 
(dashed line), Hme (dotted line), and Hmd (dot-dashed line). 
Both completions exhibit a strong spurious peak near the di- 
rection w = 0.9564, where the true spectrum has an extremely 
small local maximum. This again shows that completion 
(especially ME completion) can self-induce ambiguity where 
none actually exists; because of this phenomenon, ambiguity 
resolution is a necessary part of the DOA estimation procedure. 
The existence of a scenario-dependent bias leaves room for an 
unconventional asymptotic accuracy behavior. 

The second half of our simulation results deals with the sto- 
chastic case, where the statistical performance of our techniques 
can be established. Tables II and III summarize the results of 
stochastic experiments for the co-array-manifoldly ambiguous 
(di,W4) scenario and the unambiguous (di,u/ 4a ) scenario, 




-1 -0.75-0.5-0.25 0 0.25 0.5 0.75 1 
DOA(w) 

Fig. 3. MUSIC pseudo-spectra of H cxact , Hme, and H MD for the 
unambiguous scenario w 4a = (—0.8, -0.3, 0.1, 0.4]. 



TABLE II 

Results of Stochastic Simulations for the Fig. 1 Scenario, Each with 
1000 Monte-Carlo Trials and N Snapshots 





N = 


= 50 


N = 


200 


N = 


1000 






Hme 


Hmd 


Hme 


Hmd 


Hme 


Prob of correct identification 


0.520 


0.754 


0.502 


0.887 


0.493 


0.997 


Unrefined bias 


0.0182 


0.0017 


0.0126 


0.0007 


0.0097 


0.0003 


Unrefined standard dev 


0.0116 


0.0011 


0.0065 


0.0003 


0.0032 


0.0001 


Refined bias 


0.0086 


0.0071 


0.0028 


Refined standard dev 


0.0143 


0.0055 


0.0021 


Cramer- Rao bound 


0.0157 


0.0078 


0.0035 



TABLE III 

Results of Stochastic Simulations for the Fig. 3 Scenario, Each with 
1000 Monte-Carlo Trials and N Snapshots 





N = 


= 50 


N = 


200 


N = 


1000 




Hmd 


Hme 


Hmd 


Hme 


Hmd 


Hme 


Prob of correct identification 


0.909 


0.934 


0.994 


0.999 


1.000 


1.000 


Unrefined bias 


0.0173 


0.0260 


0.0103 


0.0220 


0.0063 


0.0200 


Unrefined standard dev 


0.0093 


0.0063 


0.0052 


0.0036 


0.0032 


0.0019 


Refined bias 


0.0080 


0.0035 


0.0015 


Refined standard dev 


0.0070 


0.0022 


0.0006 


Cramer-Rao bound 


0.0055 


0.0027 


0.0012 



respectively. Each experiment consisted of 1000 Monte Carlo 
trials. Three different sample sizes (N - 50, 200, and 1000) 
were simulated, each for ME and MD completion. The tables 
give the sample probability of correct identification (i.e., 
ambiguity resolution), the bias and standard deviation for the 
sample of worst DOA estimate in each trial, the same measures 
for the worst DOA estimate after ML refinement, and, finally, 
the CRB. 

In the case of Table II, we see that (unrefined) DOA estima- 
tion accuracy provided by ME completion is very much better 
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than that given by MD completion, but more importantly, it 
is significantly better than the corresponding CRB. Moreover, 
local ML refinement merely degrades the ME initialization 
accuracy to a level comparable with the CRB. For the unam- 
biguous scenario of Table III, we find that this remarkable 
phenomenon does not occur; instead, ME completion generates 
DOA estimates that are significantly biased and worse than the 
CRB, as we would have expected. We are forced to conclude 
that ME completion for the Table II scenario yields DOA 
estimates that are known in statistical estimation theory as 
super-efficient [18]. 

Let us now turn our attention to the probabilities of correct 
initialization (ambiguity resolution) for these two scenarios. In 
the first case, despite the extremely accurate super-efficient ME 
estimates, the probability of correctly resolving the ambiguity 
is merely quite good, increasing to 99.7% for N = 1000 snap- 
shots. This reveals that the diagonal-fitting method can be very 
sensitive to DOA estimation errors. For this reason, MD com- 
pletion (that is significantly less accurate than ME completion in 
this case) has no chance to resolve ambiguity; for all our exper- 
iments, the probability of correct identification was essentially 
equal to one half, and the power estimates for the first (true) 
and fifth (spurious) source averaged around 50 each. Only when 
the sample size was increased to N = 10 4 did this probability 
increase to 67.5%, whereas only for N = 5 x 10° were all 
1000 trials correctly identified. For the unambiguous Table III 
scenario, both completion techniques demonstrate a very high 
probability of correct identification for a very modest sample 
volume, despite significant DOA initialization inaccuracies. 

Thus, we may conclude that the NLA di and AGS wags in- 
troduced in [7] has two additional remarkable properties. First, 
DOA estimates produced by ME completion are super-efficient. 
Second, due to the above example, we may introduce the con- 
cept of unstable identifiability since neither local nonidentifi- 
ability (rank deficiency of the Fisher information matrix) nor 
pointwise nonidentifiability apply in this case. Here, we have 
demonstrated the ability to unmistakably resolve this ambiguity 
for the exact covariance lags and DOA's; however, in a small 
neighborhood of these values, the solution of the LP problem 
bifurcates, and therefore, the problem is identifiable only for an 
extremely large sample size (which is not always possible in 
practice). This phenomenon is well-known in the field of per- 
turbation analysis in mathematical programming [26]. 

Thus, any identifiable scenario that is unstable in some small 
neighborhood of the true parameters (r,0*) may be called un- 
stably identifiable since for practical applications with a finite 
sample size N, this property actually means that the scenario is 
nonidentifiable unless N is several orders of magnitude larger 
than in stably identifiable scenarios. 

It might be interesting also to comment on the feasibility sta- 
tistics of p.d. completion for the incomplete augmented sample 
covariance matrix. For the Table III scenario, for example, the 
feasibility condition H(z) > 0 did not hold in about 24% of 
our 1000 trials for N = 50 snapshots, and hence, the modi- 
fication of the specified covariance lags as described in Sec- 
tion III-C was applied. For N = 200, this probability de- 
creased to 1 1%, whereas for N = 1000, a single infeasible trial 
was found. For the more sensitive Table II scenario, the prob- 
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Fig. 4. MUSIC pseudo-spectra of H exac t, H M e, and H MD for the 
unambiguous scenario w 4b = [—0.8, —0.3, 0.1, 0.3] and "GIMP-identifiable" 
array in- 



abilities were not surprisingly higher: approximately 46% for 
N = 50, 200, and 1000, dropping to 26% only for N = 10 4 . 
Interestingly, these probabilities are quite close to those of in- 
correct identification produced by MD completion in this case. 

Our final example deals with the quasi-fully augmentable and 
GIMP-identifiable array 

d 2 = [0,1.09,3.96,5.93] , (58) 

that has the potential to unambiguously identify up to six uncor- 
rected sources. We will first demonstrate that MD completion 
can be an improvement on ME completion by considering the 
m = 4 scenario 



w Ab = [-0.8, -0.3,0.1,0.3]. 



(59) 



Of course, M = 4 is not sufficiently large to expect extreme 
consequences of the augmented covariance matrix incomplete- 
ness; however, by beginning matrix completion with the true 
covariance matrix R eX act> the MUSIC pseudo-spectrum of 
Hmd (dot-dashed line in Fig. 4) clearly demonstrates an accu- 
rate and consistent identification that is very close to the ideal 
pseudo-spectrum MUSIC (H exact ) (dashed line), whereas 
MUSIC (H me) (dotted line) has only three peaks. Obviously, 
the initialization approach we have proposed cannot work 
properly via ME completion for this scenario. This emphasizes 
that neither ME nor MD completion is an ideal technique 
for DOA estimation and that further work into completion 
philosophies is necessary. 

Finally, we consider the equispaced six-source superior sce- 
nario 

d 2 =[0,1.09,3.96,5.93] 

w 6 = [-0.80, -0.47, -0.14, 0.19, 0.52, 0.85]. (60) 

The first block of data in Table IV presents the results of sto- 
chastic simulations, each involving 1000 Monte Carlo trials for 
20 dB SNR and N = 200 snapshots. Here, we are only in- 
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TABLE IV 

RESULTS OF STOCHASTIC SIMULATIONS FOR INITIALIZATION IN THE UNAMBIGUOUS SlX-SOURCE SCENARIO AND A ' 'GIMP- IDENTIFIABLE" ARRAY FOR VARIOUS 
INTERSOURCE CORRELATION COEFFICIENTS a; 1000 MONTE-CARLO TRIALS AND 200 SNAPSHOTS. CRB FOR a = 0 IS EQUAL TO 0.0032 





a = 


= 0 


a = 


0.1 


a = 


0.5 


Of = 


1.0 










Hme 


Hmd 




Hmd 


Hme 


Unrefined bias 


0.0107 


0.0106 


0.0108 


0.0112 


0.0126 


0.0153 


0.0155 


0.0172 


Unrefined standard dev 


0.0041 


0.0016 


0.0038 


0.0015 


0.0015 


0.0005 


0.0006 


0.0003 



terested in comparing initialization accuracy between ME and 
MD completions since the asymptotic efficiency of the local 
ML refinement procedure has already been demonstrated. As 
before, the bias and standard deviation for the worst DOA esti- 
mate is given. This scenario with well-separated sources is al- 
most equally well initialized by either completion. 

The final important issue that needs to be addressed is related 
to our basic assumption that there is no intersource correlation, 
i.e., that the matrix P in (4) and (7) is purely diagonal. In any 
practical application, there might be some residual nonzero in- 
tersource correlation that should be tolerated by our approach 
without a significant degradation in performance. In order to in- 
vestigate this issue, we have repeated the stochastic experiments 
for the last six-source scenario, modifying the intersource cor- 
relation matrix in the following way: 

p=10 s^/io (/m + alm) (61) 

where l m is the m-variate unit matrix fora = 0.1, 0.5, 1.0. This 
model represents m independent sources and their fully corre- 
lated component. We can see from Table IV that the initializa- 
tion accuracy does not dramatically degrade as the intersource 
correlation increases. 

We deliberately have not performed local DOA refinement 
for both uncorrected and correlated trials (for different rea- 
sons). It is important to emphasize that augmentation relies on 
the uncorrelated-source model, whereas the refinement proce- 
dure could, in principle, incorporate any specific model for the 
intersource correlation matrix P. Since augmentation is here 
used for initialization only, it would be straightforward to in- 
troduce modifications incorporating such a model. 

Computationally, both ME- and MD-completion techniques 
are very efficient, especially compared with the computationally 
intensive one-dimensional (1-D) local ML search. 

VII. Summary and Conclusions 

This paper demonstrates that the DOA estimation problem 
for noninteger nonuniform linear arrays may be successfully 
solved for a superior number of identifiable uncorrected 
Gaussian sources. We have formulated necessary and sufficient 
conditions for both local and pointwise nonidentifiability. 
The former means that the Fisher information matrix is 
rank deficient, whereas the latter means that the m-source 
covariance matrix is identical to some similar mixture of m 
independent sources with at least one different DOA. Pointwise 
nonidentifiability is defined by the algebraic properties of the 
particular array geometry d and DOA scenario w. Since one 
of the necessary conditions for this type of nonidentifiability 
is manifold ambiguity in the co-array, the Proukakis-Manikas 



technique [7], [8] for computing ambiguity generator sets 
(AGS's) is used to list potentially nonidentiflable scenarios 
for any given array in order to test each of them against the 
pointwise nonidentifiability condition. 

Both of these types of nonidentifiability are strict in the sense 
that they occur for true values of DOA's and covariance lags. 
By an example, we have also established the existence of a dif- 
ferent type of unstable identifiability, which is important when 
considering practical applications. This occurs in scenarios that 
are theoretically identifiable, i.e. for true values, and where our 
proposed method can perfectly resolve ambiguity (again for true 
DOA's and covariance lags) but where for very small estimation 
errors, the LP solution bifurcates, meaning that the scenario is 
identifiable only for an extremely large sample size (which is 
not always possible in practice). It seems to be more than co- 
incidental that the same example of unstable identifiability also 
demonstrates super-efficient DOA estimates with one of the ma- 
trix completion methods we have discussed. Nevertheless, even 
this remarkable DOA estimation accuracy, which is far better 
than predicted by Cramer-Rao bound analysis, does not result 
in near-perfect ambiguity resolution due to the great sensitivity 
of the problem. 

These phenomena have been demonstrated using a rather un- 
usual array (noninteger but still highly redundant), introduced 
in [7] to illustrate high-order ambiguity. From this viewpoint, 
GIMP-identifiable arrays are introduced here as a preferred 
geometry; this class does not satisfy the co-array manifold 
ambiguity condition, and therefore, such arrays are globally 
identifiable (for independent Gaussian sources), modulo the 
unproven completeness of the AGS-generating technique. For 
quasi-fully augmentable (noninteger) geometries that can be 
treated as slightly perturbed versions of some corresponding 
fully augmentable (integer) array, we conjecture (for indepen- 
dent Gaussian sources) that for sufficiently small perturbations, 
such geometries are also strictly globally identifiable (regard- 
less of the AGS-generating completeness). 

Our approach to DOA estimation for noninteger linear arrays 
with an equal or greater number of uncorrelated sources than 
sensors (the superior case) essentially consists of three steps. In 
the first initialization stage, we associate our problem with some 
completion of the partially specified Hermitian covariance ma- 
trix corresponding to a particular virtual array. We have proposed 
that the co-array is the logical choice for this virtual array. The 
specified elements of the Hermitian matrix are estimated using 
the direct augmentation approach (DAA) applied to the standard 
direct data covariance (DDC) matrix. Two ad hoc completion 
criteria have been discussed: maximum-entropy (ME) and min- 
imum-deflective (MD). For ME completion, a unique optimum 
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solution is guaranteed for any feasible initial set of specified 
covariance lags. A similar uniqueness is guaranteed for MD 
completion, which was devised for slightly perturbed integer 
array (quasi-integer) geometries but may be applied in all cases. 

The feasibility condition, which is evidently satisfied by the 
existence of the true values of the specified covariance lags, gen- 
erally does not hold in the practical case of sample covariance 
lags. We used computationally efficient interior-point convex 
programming routines to modify these stochastic covariance 
lags (if necessary), resulting in the global extremum of this 
sample covariance matrix completion problem. 

These two completion methods provide a reasonable initial- 
ization accuracy in many cases; nevertheless, it is important to 
emphasize that both these ad hoc criteria are not always consis- 
tent with the maximum-likelihood (ML) principle and may be 
inappropriate for DOA estimation in some cases. For the ME cri- 
terion, this is well known since the famous ME (Burg) spectrum 
can be poorly consistent with DOA estimation. We have demon- 
strated that ME completion has, apart from a possible DOA bias, 
the potential to increase the signal subspace, i.e., to create a 
self-induced ambiguity or, on the contrary, to increase the "white 
noise level" (minimum eigenvalue) to submerge poorly sepa- 
rated sources in the new "noise level." Both these phenomena 
are direct consequences of its determinant-maximizing nature. 

The nonoptimality of both these initializing completion 
methods has made it necessary to improve the accuracy of 
the DOA estimates by introducing two subsequent refinement 
steps. The purpose of the first refinement procedure is to re- 
solve possible co-array manifold ambiguity and/or self-induced 
ambiguity. For all (stably) identifiable scenarios, the proposed 
linear programming routine is able to resolve these ambiguities 
with a reasonably high probability. 

The second refinement step is applied to the set of now un- 
ambiguous DOA's and implements a local ML search in the 
neighborhood of the initial DOA estimates. This is necessary 
to remove a possible bias introduced by completion and to en- 
sure that the final accuracy is close to the CRB. The latter is 
important because the DAA is known to be far from asymptoti- 
cally optimal [11]. Even so, the biased nature of the DOA esti- 
mates produced by ME completion has left room for super-effi- 
cient DOA initialization in one particular presented scenario. In 
less exotic examples and, specifically, for GIMP-identifiable ar- 
rays, our approach demonstrates almost CRB-optimal final ac- 
curacy. Finally, we have also seen that the augmentation tech- 
nique, which relies on the uncorrected sources model, is quite 
robust with respect to a finite amount of intersource correlation. 
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